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Abstract. We study the phase diagram and critical behavior of the two-dimensional 
square-lattice fully frustrated XY model (FFXY) and of two related models, a lattice 
discretization of the Landau-Ginzburg- Wilson Hamiltonian for the critical modes of the 
FFXY model, and a coupled Ising-XY model. We present a finite-size-scaling analysis 
of the results of high-precision Monte Carlo simulations on square lattices L x L, up 
to L = O(10 3 ). 

In the FFXY model and in the other models, when the transitions are continuous, 
there are two very close but separate transitions. There is an Ising chiral transition 
characterized by the onset of chiral long-range order while spins remain paramagnetic. 
Then, as temperature decreases, the systems undergo a Kosterlitz-Thouless spin 
transition to a phase with quasi-long-range order. 

The FFXY model and the other models in a rather large parameter region show 
a crossover behavior at the chiral and spin transitions that is universal to some 
extent. We conjecture that this universal behavior is due to a multicritical point. 
The numerical data suggest that the relevant multicritical point is a zero-temperature 
transition. A possible candidate is the 0(4) point that controls the low-temperature 
behavior of the 4- vector model. 
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1. Introduction 

In the last few decades there has been considerable interest in the nature of the 
phase diagram of the two-dimensional (2-d) fully frustrated XY (FFXY) model, whose 
Hamiltonian is [Tj 



Here J > 0, the angle variables 9 X are defined on the sites of a regular 2-d lattice, the 
summation is over nearest-neighbor pairs, and the quantities A xy are fixed and satisfy 
the constraint ^2 A X y = ft around each plaquette. The FFXY model is experimentally 
relevant for Josephson-j unction arrays in a magnetic field. It should describe the 
superconducting-to- normal transition at half a flux quantum per plaquette, see, e.g., 
Refs. El El and references therein. 

As a consequence of frustration, the ground state of the FFXY model presents an 
0(2)®Z 2 degeneracy PQ. While the 0(2) degeneracy is the usual one and is related to 
the 0(2) global invariance of the Hamiltonian that is broken in the low-temperature (LT) 
phase, the additional Z 2 degeneracy [HE]- A variable that distinguishes the different 
ground states is the chirality [I], which is defined by 



where U n is a lattice plaquette and n is the dual-lattice site at its center. At 
zero temperature there are two degenerate ground states related by a global flip of 
the chiral variables \ n . Therefore, the staggered chirality magnetization defines an 
order parameter related to the chiral Z 2 symmetry, which competes with the spin 
magnetization in determining the phase diagram and critical behavior of the FFXY 
model. Much work has been dedicated to the study of the phase diagram of the FFXY 
and related models that have the same ground-state degeneracy |H El HI IH1 EE CHI EH UH 



uni nil hhi cni nzi hhi un 1201 1211 1221 123 1231 1211 1201 inzi 12a 123 eoi eh E21 esi eh E3 ina 
cm ehi csa m3 sh S21 sni eh ssi inn sn shi sni ism eh E2i eh eh esi ism eh ehi isni Eim 

EHEHEHEHEHinnilSHinHlinnilZnilZHIZHIZSIHIZSl- In spite of that, there is not yet 



a general consensus on the critical behavior of these systems. In this paper, we address 
again this issue and, by performing simulations on very large lattices and by means of 
a careful finite-size scaling (FSS) analysis of the results, we provide what we believe is 
the definite answer on this problem. 

Two scenarios have been put forward for the critical behavior of the FFXY and 
related 2-d models. In one of them the system undergoes two transitions as temperature 
decreases. First, chiral modes undergo an Ising-type transition at T c h, characterized by 
the onset of chiral long-range order while spins remain paramagnetic. Then, at a lower 
temperature T sp , spin modes exhibit a transition to a phase with quasi-long-range order, 
which coexists with the long-range order of the chiral modes. This second transition 
is expected to belong to the same universality class as the Kosterlitz-Thouless (KT) 
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transition in the standard 2-d XY model [7B1 [77] . In this scenario spin and chiral modes 
decouple and order at different temperatures. This scenario is supported by arguments 
based on a kink-antikink unbinding picture [OBI El]- A second possibility is that spin and 
chirality order at the same temperature, T c , where both chiral long-range order and spin 
quasi-long-range order set in simultaneously. In this case spin and chiral modes may 
be coupled at the transition and may give rise to a qualitatively new critical behavior. 
Therefore, one may observe values of the chiral critical exponents that differ from the 
Ising ones, v = 1 and r] = 1/4. Field-theoretical calculations within the corresponding 
Landau- Ginzburg- Wilson (LGW) theory ED] support the existence of such a new 
universality class. 

Previous results are presented in Table ^ We report results for the fully frustrated 
XY model on the square lattice (FFXY sq ) and on the triangular lattice (FFXY tr ), 
and also for some generalizations: a model with next-to-nearest-neighbor interactions 
(FFXY nn+nn „), a nearest-neighbor model with Villain Hamiltonian, and XY models 
in which frustration is induced by the competition of nearest-neighbor and next- 
to- nearest-neighbor interactions (FXY j lt j 2 ) or by a zig-zag pattern of ferromagnetic 
and antiferromagnetic couplings (FXY ZZ ). We also include results for some related 
models that have the same ground-state degeneracy: the fractional-charge Coulomb-gas 
(CG) model that is obtained from the FFXY model by the usual duality mapping, a 
model with two coupled XY fields (2cXY) or with an Ising and XY field (IsXY), the 
antiferromagnetic XXZ model with easy-plane exchange anisotropy on the triangular 
lattice (FFXXZt r ), the frustrated antiferromagnetic six-state clock model on a triangular 
lattice (FA6SC), the 19- vertex model, the quantum ladder of Josephson junctions 
(QLJJ), the solid-on-solid model coupled to the Ising model (SOS-Is), and the 
corresponding LGW (p 4 theory. In Table^we report the method that has been employed 
to investigate the phase diagram, the number of transitions that are observed (whenever 
there are two transitions, we also give the difference of the critical temperatures) and the 
estimates of the critical exponents at the chiral transition. Most of the results have been 
obtained by Monte Carlo (MC) simulations. But we should also mention transfer-matrix 
(TM) and real-space renormalization-group (RG) calculations, as well as perturbative 
field-theoretical (FT) analyses. Note that most MC studies focussed on the finite- 
size behavior at criticality; only in a few cases jHOl EM ESI H7J was the temperature 
dependence investigated in the large-volume limit. In some cases the critical behavior 
was investigated by studying the nonequilibrium relaxation (NER) at criticality. 

The results summarized in Table ^ are rather contradictory. The most recent MC 
simulations favor the existence of two transitions. If this is correct, the results of the MC 
simulations observing only one transition may be reasonably explained by noting that 
the two transitions are very close and that high accuracy is required to disentangle them. 
If there are two transitions, chiral and spin modes decouple and thus one expects Ising 
behavior at the chiral transition. But this is not supported by numerical simulations 
that have found v ~ 0.8 in most cases. To explain this result, one may conclude that 
there is a single transition that belongs to a new universality class. In this case, the 
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observation of two very close transitions might be explained by uncontrolled systematic 
errors in the analysis of the MC data. Of course, it is also possible, as discussed in 
Refs. (101 H3 EZ1 , that the departure of the chiral exponents from the Ising values is 
due to a slow crossover towards the Ising asymptotic behavior, somehow caused by the 
interaction with the spin modes, which will eventually give rise to a KT transition. 
However, this argument does not explain why the same estimate v w 0.8 is obtained in 
several different models. Finally, there is a third, less conventional possibility: There 
are two transitions, but the chiral transition does not belong to the Ising universality 
class for some unknown reasons. 

In this paper we study the phase diagram of the 2-d square-lattice FFXY and of 
two related models, a coupled Ising-XY (IsXY) model and a 4 model obtained from a 
straighforward lattice discretization of the LGW Hamiltonian for the critical modes of 
the FFXY model, see, e.g., Refs. j2HlE]- We present MC simulations on square lattices 
L x L, up to L = O(10 3 ). The phase diagrams and critical behaviors are obtained by 
means of a FSS analysis of the MC results. Short reports already appeared in Ref. [7%j . 
In the LT phase the Z 2 chiral symmetry is broken and the spin degrees of freedom show 
the same quasi-long-range order as in the 2-d XY model. We conclusively show that the 
square-lattice FFXY model undergoes two very close but separate transitions: a KT and 
then, as temperature increases, an Ising transition with 5 = (T ch — T sp )/T ch = 0.0159(2). 
The same transition pattern is observed in the 4 and IsXY models when the transitions 
are continuous. 

Beside confirming the two-transition scenario, we have also observed an unexpected 
crossover behavior that is universal to some extent. In the FFXY model and in the 4> 4 
and IsXY models in a large parameter region, the finite-size behavior at the chiral and 
spin transitions is model independent, apart from a length rescaling. In particular, the 
universal approach to the Ising regime at the chiral transition is nonmonotonic for most 
observables, and there is a wide region in which the finite-size behavior is controlled 
by an effective exponent v e g ~ 0.8. This occurs for L < ^i c \ where £| is the spin 
correlation length at the chiral transition, which is usually large in these models; for 
example, ^ = 118(1) in the square-lattice FFXY model. This explains why many 
previous studies that considered smaller lattices always found v « 0.8. This universal 
behavior may be explained by the presence of a multicritical point, where chiral and 
spin modes are both critical. As far as its nature is concerned, our numerical data at the 
chiral transition suggest a zero-temperature multicritical point. In this possible 
candidate is the 0(4) multicritical point that is present in the 4 lattice model and 
controls the low-temperature phase of the 4-vector model. 

This paper is organized as follows. In Sec. El we define the models investigated in 
this paper. The thermodynamic quantities considered in this study are introduced in 
Sec. 01 A brief presentation of our MC simulations is given in Sec. HI Sec. El is dedicated 
to the study of the LT phase. In Sec. |H]we discuss the critical behavior of the square- 
lattice FFXY model. In Sees. and El we investigate the phase diagram of the 4 model 
and of the IsXY model respectively. In Sec. M we discuss the crossover behavior at the 
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Table 1. Overview of results for the FFXY and related models. The shorthands used 
in the second and third column are explained in the text. The fourth column reports 
the number of transitions observed, and, in the case of two transitions, the relative 
difference between their temperatures, i.e. 6 = (T crl — T sp )/T cri . The fifth column gives 
the estimates of the critical exponents associated with the chiral degrees of freedom. 
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chiral transition. Finally, in Sec. El we summarize the main results of the paper and 
draw our conclusions. In |Appendix A| we provide some details on the algorithms used 
in the MC simulations. In |Appendix B we report some results for the LT phase of the 



2-d XY model which are used in the paper. 
2. Models 

The Hamiltonian of the square-lattice FFXY model is 

'HfFXY = — J /]jxy S x ■ S y , (3) 

where the two-component spins s x satisfy s x ■ s x = 1, j xy = 1 along all horizontal 
lines, while along vertical lines ferromagnetic j xy = 1 and antiferromagnetic j xy = — 1 
couplings alternate. Here (as in the following models) J plays the role of inverse 
temperature, the Gibbs probability being simply proportional to e~^ FFXY . 
We also consider the 4 model on a square lattice. The Hamiltonian is 

(xy),i i,x 

+ 2(U + D)J2<l>lJlx, (4) 

X 

where % = 1,2, (f>i jX is a real two-component variable, the first sum goes over all nearest- 
neighbor pairs, and 0f = 4>i-4>i- Hamiltonian TL^ describes two identical 0(2)-symmetric 
models coupled by an energy-energy term. This model is a straightforward lattice 
discretization of the LGW Hamiltonian 

^lgw= / d d x{\ £ [(dM' + r^ +^u (J2 ^) 2 + i^}(5) 

J o=l,2 ' o=l,2 

(0 a , a = 1, 2, is a two-component vector), which can be obtained by applying a Hubbard- 
Stratonovich transformation to the FFXY model, dropping terms of order higher than 
four |2"§1 1§1 Therefore, the lattice 4 model (jlj) represents an effective ferromagnetic 
theory that is expected to describe the critical modes of the FFXY model. This approach 
was already used to investigate the critical behavior of antiferromagnets on a stacked 
triangular lattice [T9] . 

It is interesting to note that the symmetry of Hamiltonian (jlj) is larger than that 
of the FFXY model. Indeed, the Hamiltonian is symmetric under separate rotations 
of the two fields and under the Z2 transformation <f>\ <-> 02- Therefore, the overall 
symmetry is [0(2)©0(2)]®Z 2 . For D > the ground state corresponds to 4>\ = and 
02 7^ or the opposite, and thus it has the same degeneracy as the ground state of 
the FFXY model. Indeed, the lowest-energy configurations are determined once one 
fixes which field does not vanish — thereby breaking the Z 2 interchange symmetry — and 
the direction of the nonvanishing one. Thus, the relation between the FFXY model 
and the 4 Hamiltonian is not fixed by the symmetry of the original Hamiltonian but 
rather by the ground-state degeneracy group, which is the quotient of the symmetry 
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groups of the model and of the ground state. This is consistent with the paradigm that 
relates the universality class to the symmetry breaking pattern. Note that the extended 
symmetry of the LGW Hamiltonian is related to the truncation of the Hamiltonian to 
fourth order in the fields and it is lost if higher-order terms are added. More precisely, 
the Hubbard- Stratonovich calculation gives rise to <ft 6 and <ft 8 terms on the triangular 
and square lattice respectively, which break the 0(2) ©0(2) symmetry down to 0(2) 
These terms are irrelevant close to four dimensions but it is far from clear that the 
same holds in two dimensions. In any case, they do not change the relevant symmetry 
breaking pattern. 

Due to the ground-state structure, for J — > oo the field-interchange symmetry 
0i ^ 02 is broken and thus this symmetry is the analog of the Z 2 chiral symmetry of 
the FFXY model. The corresponding order parameter is 

C x = 01, - 4>l x - (6) 

For D = model (jlj) is 0(4) symmetric. Therefore, it does not have any transition at 
finite temperature. Criticality is observed only for J —>■ oo. In this limit the correlation 
length £ increases exponentially, i.e. £ ~ e cJ ; see, e.g., Ref. 
We also consider the Ising-XY (IsXY) model 

J 



^ (1 + a x a y ) s x ■ s y + Ca x a y 



(7) 



^IsXY — — / ' J 

where a x = ±1 and the two-component spins s x satisfy s x ■ s x = 1. Here s x and o~ x 
correspond to spin and chiral variables, respectively. Note that, by performing the limit 
U — > oo and then D — > oo in the 4 model one recovers the IsXY model for C — 0. 
In this case, the variables a x and s x are related with those of the 4 model by 

01,* = |(1 +0- x )s x , 2 ,a; = |(1 - (T X )S X . (8) 

Models with C ^ can also be recovered from a 4 Hamiltonian. It is enough to add 
an energy-energy nearest-neighbor hopping term Yllxy) (0i,a;0l,3/ + 02,x02,y)- 

Apparently, the IsXY model is only invariant under the group 0(2)®Z 2 . However, 
its relation with the 4 model indicates that the symmetry is larger. Indeed, for any 
value of C, the model is invariant under the 0(2)©0(2) nonlinear transformations 

si = [|(1 + <t x )RW + |(1 - o- x )R^}s x , 

< = <Tx, (9) 

where and R^ are 0(2) rotation matrices. 

It is possible to argue that in the IsXY model chiral and spin modes cannot be 
critical at the same value of J and C, see, e.g., Ref. |[)8] . Let us assume that, for fixed 
C, the model undergoes a continuous phase transition at J = J c h, where the chiral 
correlation length is infinite, i.e. where a correlations are critical. Let us consider first 
the dynamics of the spins s x for a fixed Ising {c x } configuration. Spins a with the 
same sign form geometrical clusters and spins s x belonging to different geometrical (not 
Fortuin-Kasteleyn) clusters do not interact (either directly or indirectly) because of the 
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prefactor (1 + o~ x a y ) in the hopping term of the spin variables. Thus, the behavior of 
the spins s x is completely determined by their behavior on each single cluster. It is easy 
to convince oneself that J = J c h is the percolation threshold of the clusters. Indeed, 
if the clusters percolate at J* < J c h, for J = J* Ising spins would already be critical, 
since in two dimensions percolation of the geometrical clusters implies Ising criticality. 
Now, at the percolation threshold we expect the clusters to have a fractal dimension 
that is smaller than 2. This is true for the Ising transition [EHIE2]; we assume here that 
it holds in general. Thus, fields s x interact as on a system with d < 2 and thus cannot 
be critical. Therefore, at the chiral transition the spin correlation length must be finite, 
and, moreover, J c h < J sp . As we shall see, the MC results that we shall present provide 
support to this argument. 



3. Definitions and notations 

3.1. The fully frustrated XY model 

In the square-lattice FFXY model the ground state is only invariant under translations 
of two lattice spacings. Therefore, we divide the lattice into four sublattices, so that 
the four sites of each plaquette belong to different sublattices. Then, we define the 
spin correlation function G s {x) = (sq ■ s x ) only for x = (2n,2m), n,m integers, and its 
Fourier transform as G s (q) = Yl x e tqx G s (x) where the sum goes over x = (2n, 2m). The 
corresponding susceptibility Xs an d second-moment correlation length £<, are given by 

x=(2n,2m) 

1 G a (0)-G a (q) 
^ s 4sin 2 (g min /2) G s (q) ' 

where q = (g m in,0), and g min = 2n/(L/2). We shall consider two RG invariant ratios 
related to the continuous spin modes: 

L r> — 

where fx is the magnetization corresponding to one of the four sublattices, i.e. 

x=(2n,2m) 

V = L 2 is the volume, and fj, 2 = ft - jl. We also define the helicity modulus T. For this 
purpose we introduce a twisted term in the Hamiltonian. More precisely, we consider 
the nearest-neighbor sites (x, y) with x\ = L, yi = 1, and x^ = y2, and replace the term 
s x ■ s y in Hamiltonian (JHJ) with 

s x ■ R^s y = 4 X) (4 1} cos V 5 + 4 2) sin V 9 ) + s x ] ( s y 2) cos V 3 - s x 1] sin cp) , (14) 



R, ^ B 3 = (12) 
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T 



(15) 

tp=0 



where R^ is a rotation of an angle tp. The helicity modulus is defined as the second 
derivative of the free energy with respect to <p at p = 0: 

d 2 In Z{<p) 
dip 2 

Chiral modes are related to the Z2 symmetry that is broken by the ground state. As 
discussed in the introduction, a good order parameter is (we have dropped the sign here 
and thus C n is no longer a "spin" susceptibility; this is of course irrelevant) 

C n = jxysm(9 x - 9 y ), (16) 

where n is the dual lattice site at the center of Il n . We consider the staggered correlation 
function 

G c (n) = (-ir^(C C n ) c , (17) 

the chiral susceptibility Xc, and the second-moment correlation length £ c . They are 
defined as in Eqs. (110)1 and (JUJ, though in this case sums are extended over the whole 
dual lattice and g m i n = 2n/L. Analogously, we define the RG invariant quantities 

Rc = |, B c EE (18) 

where the staggered magnetization is defined as 

^ = v X (-ir +n2 ^. (19) 

n=(ni,n2) 

3.2. The (j) 4 model 

In the </> 4 model we define the hopping energy density and specific heat as 

E=^(H h h C ee i ((Hi) - (H h ) 2 ) , (20) 

where 

{xy) 

The spin two-point correlation function G s (x) is defined as 

G s (x) EE (0 1>o ■ M + 02,0 ■ 4>2,x)- (22) 

The corresponding susceptibility Xs and second-moment correlation length £ s are given 
by 

x< s E G , w , g s i (23) 

V ' 4 S in 2 (,, mil ,/2) G ,(q) 

where G s (q) is the Fourier transform of G s (x), q = (q m in, 0), and q m \ n = 2n / L. We shall 
consider two RG invariant ratios related to the continuous spin modes: 
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where 

fk = y ^ A = ^i- lk- ( 25 ) 

X 

As before, we consider the helicity modulus Tj associated with each variable 0^ = 
(<j>fxi$ix)i defined analogously to Eqs. (fTljl and (fT3j) . The total helicity modulus is 
given by T = T x + T 2 . 

The "chiral" modes related to the Z 2 field-interchange symmetry are associated 
with the quadratic operator defined in Eq. (j5J. We consider the corresponding connected 
two-point correlation function 

G c (x) = (C C x ) c , (26) 

the chiral susceptibility Xc the second-moment correlation length £ c defined as in 
Eq. (j2Bj) , and the ratios 

Cc D - (/4) a^Ec. ( 27 ) 



R.r = ~r, B, 



V (/x c 2 ) 2 ' ^ v 

3.3. The Ismg-XY model 

In the case of the IsXY model (J7J) we consider the correlation functions G s (x) = (s ■ s x ) 
and G c (x) = (ero • <j x ) and the corresponding observables that are defined as in the 4 
model. For C = it is easy to verify that G s (x) and G c (x) exactly correspond to the 
correlation functions in the 4 model. Indeed, using mapping (JHJ), C x = a x , so that 
chiral correlations correspond to correlations of the Ising variables. On the other hand, 
mapping (jHJ) gives 

G s (x) = + a x o y )s x ■ s y ). (28) 

However, since Hamiltonian (jjj) is invariant under the transformations s x —>■ o~ x s x , we 
have 

(a x a y s x ■ s y ) = (s x ■ s y ). (29) 

Therefore, 

G.{x) = (s x ■ s y ). (30) 

Note that the Binder parameter Bg^, defined as in Eq. using the mapping (jHJ), is 
not the natural one in terms of the spin variables s x . Indeed, mapping (jHJ gives 

R l ((^) 2 + 2/i^ + (/iL) 2 ) 

B « ~ 4 > ( 31 ) 



where 



The conventional Binder parameter is instead^ 

B s = ((^) 2 >/(^> 2 - (33) 
In the HT phase we have B S( p = 3/2 and B s = 2 in the thermodynamic limit. 
| In terms of the original 4> variables, we have B s = B s( f, + ^(uflJ^) / Q^iit 1 !)) 2 ■ Note that B s > B^. 
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4. Monte Carlo simulations 



We perfomed MC simulations of the FFXY model (JTj), of the 4 model (JH) for U = 1 
and D in the range 1/20 < D < 99, and of the IsXY model for -8 < C < 0.3. We 
considered L x L square lattices with periodic boundary conditions, up to L = O(10 3 ). 
The phase diagram and critical behavior is investigated by analyzing the FSS behavior 
of several quantities. The total CPU time used in this study was approximately 10 CPU 
years of a single 64-bit Opteron 246 (2Ghz) processor. 

In the MC simulations of the FFXY model, we used a local algorithm based 
on a mixture of Metropolis and overrelaxed (microcanonical) updates, as suggested 
in Ref. [60J. In the case of the IsXY and 4 model the updating algorithm was 
based on mixtures of Metropolis, overrelaxed (microcanonical), and single-cluster [SHI 
updates. More details on the MC algorithms are reported in |Appendix A At the 
chiral transition, single-cluster updates have little influence on the autocorrelation times, 
while the overrelaxation updates are only able to speed up the simulation but do not 
change the dynamic critical exponent. This is always approximately equal to 2, as 
it is expected in the case of a purely local dynamics. This means that very long 
simulations are needed to obtain reliable statistics. For instance, in order to obtain 
approximately 2500 independent configurations for the FFXY model on a lattice with 
L = 1000, we used approximately 800 days of a single 64-bit Opteron 246 processor. 
Analogously, the runs on the largest lattices for the IsXY model (C = 0, L — 360, 4000 
independent configurations) and for the 4 model (D = 1/2, L = 1200, 150 independent 
configurations) took approximately 100 days in both cases. A substantial reduction of 
the critical slowing down is achieved at the spin transition and in the LT phase, because 
the overrelaxed and the cluster algorithm are very effective in dealing with the spin 
modes. 



5. The low-temperature phase 

In the high-temperature (HT) phase, in which symmetry is not broken, the spin and 
chiral correlation functions G s (x) and G c (x) decay exponentially at large distances. 
In the LT phase (J large) instead, the spin correlation function G s {x) is expected to 
decay with a power law, giving rise to quasi-long-range order. Indeed, according to the 
Mermin- Wagner theorem J84J, any magnetization breaking a continuous symmetry is 
forbidden in two dimensions. On the other hand, the discrete Z 2 symmetry may be 
broken with a nonvanishing chiral magnetization. In this section we study the main 
features of the LT phase. We shall see that Z 2 symmetry is broken and that 0(2) 
quasi-long-range order is realized in the whole LT phase. The critical behavior of the 
spin variables is controlled by a line of Gaussian fixed points that are exactly those that 
control the LT phase of the standard XY model. 

Direct evidence for the breaking of the Z 2 chiral symmetry in the LT phase is 
provided by the finite-size behavior of the chiral Binder parameter B c . In the LT phase 
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Table 2. Results for B c , R s , T, \s, and r) e g in the LT phase. In the case of the 4 
model we report (j\ a /L, Ti a , and Xla that converge to i? s , T, and Xs as L — > oo. The 
exponent ij e ff corresponds to w(L,2L) defined in Eq. I|34|). 
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</> 4 , -D = 1/2 


1.48 
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1.0922(13) 
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634.1(9) 
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1.0382(9) 
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7824(7) 
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512 
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we should have B c — > 1 + 0(1/^) when L — ► oo. MC simulations (see Tabled show 
this behavior in all models that we have considered. 

Let us discuss the 4 model in more detail. We distinguish the fields <pi X from the 
size of the corresponding quantity Q. t = J^x^fx- For each configuration, we define the 
large (small) field 0i a;X (<f) sm , x ) to be the field 0^ that has the maximum (minimum) 
value of Qi. Then, we define the corresponding two-point functions G\ a (x) = (0i a ,o • 4>ia,x) 
and G sm (x) = (0 S m,o • 4>sm,x), the magnetic susceptibilities xia and Xsm, and the second- 
moment correlation lengths and £ sm , defined according to Eq. (j23j). Our numerical 
results show that, as L increases, small and large fields effectively decouple. Indeed, Xsm 
converges to a constant, while xia diverges as L 2 " 71 . Therefore, in the large-L limit, since 
Xs = Xia. + Xsm, we have X\JXs = 1 + 0(L~ 2+V ). Fig.Qshows estimates of Xs, Xsm, and 
Xi a for U = 1, D = 1/2, and J = 1.48 in the LT phase [as we shall see, for D = 1/2 there 
are two very close transitions, an Ising transition at J c ^ = 1.4668(1) and a KT transition 
at J sp = 1.4704(2)]. It is evident that the large field is critical for any J, unlike the 
small one. We also checked the behavior of the correlation lengths: while £i a diverges 
as L, the small component £ sm remains finite. As a consequence £ia/£s = 1 + 0(L~ 2+V ), 
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Figure 1. Magnetic susceptibilities Xs, Xla, and Xam for U = 1, D = 1/2, and 
J = 1.48. The data for Xs an d xia can hardly be distinguished on the scale of the 
figure since differences are very small. 
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Figure 2. Estimates of R s vs r/ in the LT phase. The continuous line represents the 
prediction obtained by assuming that the LT phase is described by the same Gaussian 
line of fixed points as the XY model; see |Appcndix B| 



as it was the case for the magnetic susceptibility. Finally, we computed the helicity 
moduli Ti a and T sm . The helicity modulus of the small field vanishes within error bars 
starting from relatively small lattice sizes (for example, for L > 64 for J = 1.48 and 
D = 1/2). Therefore, also for this quantity we have T = T\ + T 2 « Ti a . These results 
show convincingly that the large- distance behavior in the LT phase of the 4 model 
is effectively determined by only one of the two fields 0j, with a complete decoupling 
of the other one. This reflects the fact that the Z 2 chiral symmetry is broken with a 
spontaneous magnetization Mc ^ 0. 
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Table 3. Estimates of R s = £ S /L, T, and r\ in the LT phase. The column Z/ max 
reports the size of the largest lattice used. The estimate of 77 that appears in the fifth 
column is obtained by using the finite-size behavior of Xs ~ L 2 ^ 71 . Those given in the 
last two columns are obtained by using the estimates of R s and T and the relations 
among these variables that are obtained in the XY model. 
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We now prove that the long-distance behavior of the spin correlation function is 
analogous to the one of the standard XY model, i.e. it is controlled by the same line 
of Gaussian fixed points with varying critical exponent rj. For this purpose, we show 
that, in the large- L limit, the exponent n (computed from Xs ~ L 2 ~ v ), R s = £ S /L, and 
T satisfy the same universal FSS relations as the corresponding quantities in the XY 
model. These relations, that are reported in |Appendix B[ are universal once the shape 
of the lattice and the boundary conditions are fixed: in our case we consider a L x L 
square lattice with periodic boundary conditions. 

In order to verify these relations, for each L we compute R s , T, and an effective 
estimate of r\. Since the susceptibility behaves as Xs ~ L 2 ~ v , we define 

»j(Li, L 2 ) = 2 - ln[ X (L 1 )/x(L 2 )]/ln( J L 1 /L 2 ). (34) 

In the 4 and IsXY models we use the results for %i a , £i a /£, and Ti a because they have 
smaller scaling corrections with respect to the analogous quantities \i £s/L, and T. 
Some results are reported in Table El In most of the cases, as L increases, all quantities 
approach a constant within error bars, indicating that we have reached the infinite- 
volume limit within errors. When convergence is not clearly observed within errors, we 
extrapolate the data assuming the scaling corrections that are expected in the LT phase 
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of the XY model: we perform an extrapolation of the form a + bL with ( = I/77 — 4 
for T and ( = Min(2 — 77, 1/77 — 4) for \s and R s , see Appendix B 



Infinite- volume estimates of 77 are reported in Table El Beside the values obtained 
from the magnetic susceptibility, we also report the estimates obtained from R s and T 
using the XY relations. In all cases the estimates are in good agreement. In Fig. El 
we plot ^ S /L versus 77 as obtained from the susceptibility, for several models in the LT 
phase. The results agree quite precisely with the curve obtained in the standard XY 
model, see | Appendix B Thus, we conclude that the LT phase of the FFXY, 4 , and 
IsXY models is controlled by the same line of Gaussian fixed points that are relevant 
for the XY model. 

In the XY model the LT phase becomes unstable against dissociation of vortices 
when r] = 1/4, and [EH! U L = r kt = 0.750691..., T = T KT = 0.636508.... At 
these values, a KT transition occurs and the system becomes disordered. We expect 
that a similar mechanism is realized in the models considered here, as long as the 
chiral Z2 symmetry is broken. With decreasing J, 77 increases. If 77 = 1/4 is reached, 
vortices disorder the system and thus we expect a phase transition. If the transition is 
continuous and the chiral Z2 magnetization does not vanish at the transition, then the 
phase transition is expected to be of KT type. An Ising or first-order transition should 
then occur at a smaller value of J, where also chiral variables disorder. Otherwise, if 
the Z2 magnetization vanishes, the transition should belong to a new universality class. 
Note that this possibility is excluded in the IsXY model by the argument reported in 
Sec. El If a transition occurs for 77 < 1/4, it cannot be a KT one. Thus, it must be of 
first order or a continuous transition belonging to a new universality class. In the latter 
case, however, we expect a single transition, since a new type of transition can only be 
due to the interaction of spin and chiral modes that become critical at the same value 
of J. 



6. Phase transitions in the fully frustrated XY model 

In this section we investigate the phase diagram of the square-lattice FFXY model. We 
show that there are two very close but separate transitions, a KT transition at J = J sp 
involving the spin degrees of freedom and an Ising transition at J ch < J sp associated 
with the breaking of the chiral Z2 symmetry. 



6.1. Finite-size scaling at the chiral transition 

We perform a FSS analysis using the method proposed in Ref. JHEj and also discussed in 
Ref . jHZl ■ Instead of computing the various quantities at fixed Hamiltonian parameters, 
we compute them at a fixed value of a chiral — this guarantees that we are observing the 
chiral transition — RG invariant quantity; in our specific case we choose R c = £ C /L. This 
means that, given a MC sample generated at J = J run , we determine, using standard 
reweighting techniques, the value Jr c such that R C (J = Jr c ) = R C) sx- All interesting 
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Table 4. Results at fixed R c = R\ s for the FFXY model. Estimates of the 
pseudocritical coupling Jr c , the chiral magnetic susceptibility Xc, the chiral Binder 
parameter B c , the derivatives of R c and B c with respect to J, the spin correlation 
length £ s , the derivative of R s with respect to J, and the helicity modulus T. 
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2.1642(16) 


28.38(3) 


1.1913(7) 


2.099(10) 


0.903(8) 


5.430(8) 


0.805(5) 


0.8628(14) 


12 


2.1735(11) 


56.99(7) 


1.1820(6) 


3.14(2) 


1.323(12) 


8.051(12) 


1.179(8) 


0.7833(15) 


16 


2.1797(8) 


93.09(11) 


1.1839(7) 


4.34(2) 


1.885(2) 


10.533(14) 


1.653(12) 


0.730(2) 


24 


2.1849(6) 


185.7(3) 


1.1854(8) 


6.98(4) 


3.06(3) 


15.26(2) 


2.53(2) 


0.649(2) 


32 


2.1887(3) 


303.4(3) 


1.1872(6) 


9.79(4) 


4.39(3) 


19.72(2) 


3.45(2) 


0.591(2) 


48 


2.1931(4) 
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1.1934(8) 


16.08(10) 


7.43(8) 


28.00(4) 
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0.514(2) 


G4 


2.1968(3) 


1004(2) 


1.1952(10) 


23.1(2) 


10.77(14) 


35.75(7) 


7.48(7) 


0.457(3) 
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2.2004(2) 


2038(4) 


1.2009(10) 


38.7(3) 


18.4(2) 


49.75(9) 


11.61(8) 


0.384(3) 
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2.20242(8) 


3385(4) 


1.2007(6) 


54.8(3) 


25.6(2) 


61.92(8) 


15.67(8) 


0.325(2) 


192 


2.20451(8) 


6932(14) 


1.2002(9) 


88.6(6) 


40.2(5) 


81.9(2) 


22.8(2) 


0.235(4) 
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2.20529(4) 


11469(20) 


1.1962(8) 


120.5(9) 


53.0(6) 


96.0(3) 


27.8(3) 


0.173(4) 
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72.8(1.3) 
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-0.005(6) 



observables are then measured at J = J# c ; their errors at fixed R c are determined by 
a standard jackknife analysis. The pseudocritical temperature Jr c converges to J c h as 
L — > oo. This method has the advantage that it does not require a precise knowledge 
of the critical value J c h (an estimate is only needed to fix J run that should be close to 
J c h). Moreover, for some observables the statistical errors at fixed R c are smaller than 
those at fixed J (close to Jch). This is due to cross correlations and to a reduction of the 
effective autocorrelation times. The ratio E/EfaR c of the errors at Jr c computed at fixed 
J and at fixed R c turns out to be roughly independent of L: E/E^ Rc ps 1.6, 4.3, 1.9, 1.0 
respectively for Xc B c , dB c /dJ and T. 

The value -R Ci fbc can be specified at will as long as it positive, § though R c fi x = R c (R c 
is the critical-point value) improves the convergence for L — > oo jHEllEZI- Since in the 
two-transition scenario the chiral transition is predicted to be in the Ising universality 
class, we choose R c fi x = Ri s where Ri s = 0.9050488292(4) is the universal value of £/L 
at the critical point in the 2-d Ising universality class 88J. This choice does not bias our 
analysis in favor of the Ising nature of the chiral transition. For any chosen value (as 
long as it is positive) and whatever the universality class of the chiral transition is (it 
may also coincide with the spin transition), Jr c (L) converges to J c h for L — > oo. Indeed, 
in the FSS limit R c = /r[J 1 ^(J — J c h)], where v and fn(x) depend on the universality 
class of the transition. Therefore, fixing R c is equivalent to fixing X = L x ' v \Jr c {V)— Jch]. 
Critical exponents determined by using FSS do not depend on the chosen value R c fi x , 
while other finite-size quantities do. For instance, the Binder cumulant satisfies the 

§ The method can be used with any RG invariant quantity. For instance, one could also use the Binder 
parameter. The fixed value can be taken arbitrarily between the corresponding HT and LT values. 
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Figure 3. R s = £ S /L (above) and £ s (below) at fixed R c versus 1/L and L respectively, 
for the FFXY model, the cf) 4 model at D = 1/2, and the IsXY model at C = 0. 



analogous relation B c = f B [L l / v {J - J ch )], so that B c converges to /b(X) in the limit 
L — > oo at fixed R c . By fixing R c to the critical Ising value, we will be able to perform 
an additional consistency check. If the chiral transition belongs to the Ising universality 
class, then X = (apart from scaling corrections) and we should find that any RG- 
invariant quantity converges to its critical-point value in the Ising model. 

If the transition is unique also spin variables should have a finite nontrivial limit: 
for instance, R s at fixed R c should converge to a nonvanishing constant. On the other 
hand, if the chiral transition occurs in a region in which the spin variables are still 
disordered, the spin correlation length should be finite and thus R s at fixed R c should 
converge to zero as 1/L. 

We performed simulations on quite large lattices, up to L = 1000. Results are 
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Figure 4. AB C = B c - B Is at fixed R c = R ls vs L~ 7 / 4 , for the FFXY model, the 
</> 4 model at D = 1/2, and the IsXY model at C = 0. B Is = 1.167923(5) is the value 
of the Binder parameter at the critical point in the Ising model The dotted lines 
correspond to a linear fit of the data for the largest available lattices. 



reported in Table EJ In Fig. El we show the spin correlation length £ s as a function of L. 
It clearly converges to a constant as L — > oo. We can thus estimate the infinite- volume 
spin correlation length £s at the chiral transition: ^ = 118(1). Moreover, the helicity 
modulus T at the chiral transition appears to vanish in the large- L limit, see Table EJ 
consistently with the fact that spin variables are disordered. These results definitely 
show that spin modes are not critical at the chiral transition, and therefore that, for 
J > J eh, there is a phase in which the chiral Z 2 symmetry is broken while spins remain 
paramagnetic. 

If spin and chiral modes decouple, the chiral transition should belong to the Ising 
universality class. The best evidence for that is provided by the finite-size behavior 
of the chiral Binder cumulant B c . In the case of an Ising transition, the asymptotic 
behavior of B c is expected to be [HB1 EH1 EH 

B c = B ls + 0{L-») + 0(L- 2+?? ), (35) 

where JSHJ B\ s = 1.167923(5). The first correction is due to the leading irrelevant 
operator, the second one to the analytic background terms. Since for unitary Ising 
models JESHSD! oj = 2 and 2 — 77 = 7/4, the dominant scaling corrections are those related 
to the analytic background. The estimates of B c reported in Table H] show a peculiar 

(c) 

intermediate nonmonotonic crossover: for L ^ w 118, B c is approximately constant 
and equal to 1.19-1.20, while the asymptotic regime given by Eq. ()35|) is observed only 
when L > ^ c) . In Fig. H we plot AB C = B c - B ls versus L ' . The data clearly 
approach the Ising value as L increases. A linear fit to a + bL^ 7 ^ of the data with 
L > 384 gives B c = 1.1678(16). We also mention that a linear fit of AB C to &L~ 7 / 4 for 
L > 384 gives b = 585(32) with x 2 /d.o.f. ~ 0.2, where d.o.f. is the number of degrees 



Multicritical behavior in the fully frustrated XY model and related systems 

1.5r 



19 



1/V 



eff 
0.5 



0.0 



I I I I 






A ' • 




* 




• from R * 




c 

- . from B 

c 




- a from R 




s 



100 



200 



300 



400 



500 



Figure 5. Effective exponent 1/vefi as obtained from the derivatives of R c , B c , and 
R s , for the FFXY model. The dotted line corresponds to the Ising value v = 1. 
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Figure 6. Effective exponent T] e s as obtained from Xa for the FFXY model. The 
dotted line corresponds to the Ising value r\— 1/4. 



of freedom of the fit. These results provide a rather strong evidence of the Ising nature 
of the critical behavior. 

In order to determine the critical exponents, we consider the derivatives of R c and 
B c with respect to J at fixed R c . In the large-L limit, they behave as 
dRc dB c x , 

-dT~lu~ L ■ ' (36) 

where v is the chiral exponent associated with the correlation length. Then, given 
an interval L\ < L < L2, with L2/L1 ~ 2, we compute an effective chiral exponent 
l/f e ff(Li) by fitting the available data for dR c /dJ and dB c /dJ in the interval. In some 
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cases, only two results are available and we use 

\ndS/dJ\ L2 -lndS/dJ\ Ll 

1/Ucs{Ll) = HWZo ' ( } 

where S is either R c or B c . The results are shown in Fig. |SJ The exponents l/v e g(L) 
obtained from R c and B c show a similar behavior. They are nonmonotonic and, for 
20 < L < 70, they are approximately constant and equal to 1.3, which corresponds 
to z/ e ff ~ 0.8. This is consistent with the result v « 0.8 obtained in previous studies 
that considered lattice sizes in this range [23 HI! 123 EE EH For L > = 118, the 
effective exponent l/v e s(L) decreases rapidly, undershoots the Ising value v = 1, which 
is then approached from below. In Fig. 03 we also report l/z/ e ff(L) computed by using 
R s . For small values of L, we find v e s{L) « 0.9 — this explains why spin modes have 

(c) 

sometimes been thought to be critical at the Ising transition — while for L > £g it starts 
to decrease towards the asymptotic value —1, corresponding to a finite Q at the chiral 
transition. In Fig. |U] we show the effective chiral exponent t] e s defined analogously to 
z/ e fj in terms of Xc which should behave as L 2 ~ v . The exponent is first larger than the 
Ising value, then becomes smaller, and only for L > 500 appears to approach the Ising 
value r] = 1/4. 

The critical value J c h of the hopping coupling at the chiral transition can be 
determined from the finite-size behavior of Jr c . Assuming that the transition is of 
Ising type, the large-volume behavior of Jr c is given by 

Jr c = Jch + OiL- 1 ^) + 0(L- 1 ^- 2 ^), (38) 

where the first correction is due to the leading irrelevant operators, while the second 
correction is due to the analytic background jHHl ED] • Note that in Eq. (|3~H|) we use the 
fact that the fixed value R c fi x is the critical Ising value Ri s . If R c fi x ^ R\ s , the leading 
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Table 5. Results for J sp in the FFXY model. 



Li 


Li 


from R s 


from T 


128 


256 


2.236(5) 


2.2405(12) 


256 


512 


2.2400(10) 


2.2415(5) 


512 


1024 


2.2409(3) 


2.2415(5) 



corrections are of order L~ x l v jHSUHZl- According to Eq. flUD, since to — 2 and rj = 1/4, 
the leading corrections to the infinite-volume limit are of order L _11//4 . In Fig.[7|we plot 
Jji c versus L _11//4 . In analogy with what is observed for the other chiral quantities, the 
asymptotic behavior sets in only for L ^> £g. A linear fit of the data with L > 384 
gives J c h = 2.20632(5). This estimate is significantly more precise than earlier results: 
J ch = 2.2051(24) (Ref. J ch = 2.2002(10) (Ref. [EH!), J c h = 2.212(5) (Ref. HO]), and 

J ch = 2.203(10) (Ref. [HE]). 

In conclusion, a careful FSS analysis on lattice sizes up to L = 1000 is consistent 
with an asymptotic critical behavior belonging to the Ising universality class, as expected 
within the two-transition scenario. However, before observing the asymptotic Ising 
behavior, we find a peculiar nonmonotonic crossover behavior. This is related to the 
fact that the spin correlation length £s is large at the chiral transition. For L < £^ 
spin modes appear as being critical, as also indicated by the effective exponent z/ e ff that 
is obtained from dR s /dJ, which, for L < £g , behaves as the exponent obtained from 
chiral observables, see Fig.EJ We will further discuss these crossover behaviors in Sec. El 



6.2. The spin transition 

Since spins are paramagnetic at the chiral transition, the quasi-long-range order that 
characterizes the LT phase requires a spin transition at a larger value of J. Since the LT 
phase is controlled by the same line of fixed points that are relevant for the XY model, 
and chiral modes do not play any role at the transition, it is natural to conjecture that 
it belongs to the KT universality class. 

In order to determine the critical hopping parameter J sp of the spin transition, we 
exploit the predicted asymptotic behavior of the helicity modulus T and of the ratio 
R = £/L at the KT transition. On a square lattice L x L at the critical temperature 
we have for L — > oo [85J 

0.318899454... , , 

T = 0.63650817819... H — + 39 

In L + C T v ' 

0.212430... . x 

# = 0.7506912... + -— — + 40 

\nL + C R y ' 

where Cr and CS~ are nonuniversal constants and the neglected corrections are of order 

In In Lj In 2 L. 

The critical point J sp is obtained by determining the value of J where the behavior 
of T and R s is given by Eqs. (|3*9~|) and (|4Hjl . In practice, we first obtain an approximate 
estimate of J sp by looking for the value of J at which T and R s are both close to the 
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Figure 8. Sketch of the phase diagram of the 4> 4 model for U = 1 and D > in 
the plane J~ 1 -D. The continuous, dashed, and thick continuous lines represent Ising, 
KT, and first-order transition lines. The distance between the Ising and KT lines is 
amplified; otherwise, the two transitions could not be distinguished on the scale of 
the figure. The phase diagram within the circled region is unknown. In the figure we 
show one possibility, in which the two continuous transitions connect to the unique 
first-order line. 



theoretical large-L value. Then, we perform simulations at two values of J, J\ and 3%, 
that are close to the estimated J sp , and determine T and R s for J x < J < J 2 by linear 
interpolation. Finally, for each pair of sizes L\ and L2, we require T and R s at J = J sp 
to be exactly given by the previous expressions. In this way we obtain two nonlinear 
equations for two free parameters: the critical value J sp and the constant C. Table 03 
shows the results for several pairs of lattices with L2/L1 = 2. In spite of the expected 
slow convergence, the results obtained from R s and T are consistent. This provides a 
nontrivial check of the KT nature of the transition. The helicity modulus T shows the 
least L dependence and provides our final estimate: J sp = 2.2415(5). Previous estimates 
of J sp are: J sp = 2.227(5) (Ref. [7T]), J sp = 2.2422(5) (Ref. 00]), and J sp = 2.273(10) 
(Ref. [SB])- We also estimate the infinite- volume chiral second-moment correlation length 
^ c s) at the spin transition: £c = 8.0(5). The Ising and KT transitions are very close: 
we have 

5= Jsp ~ Jch =0.0159(2). (41) 

Jch 

7. Phase transitions in the <p A model 

In this section we study the phase diagram of the 4 model for U = 1. We only consider 
the case D > 0, since only in this case is the ground-state degeneracy the same as in 
the FFXY model. We show that there is a unique first-order transition for D > D*, 
34 < D* < 49. On the other hand, for sufficiently small D, we provide evidence of two 
very close transitions, a KT transition at J = J sp associated with the breaking of the 
continuous degrees of freedom and an Ising transition at J c h < J sp associated with the 
breaking of the interchange (chiral) Z2 symmetry. In Fig. IHlwe show the phase diagram 





Figure 9. Distributions of the hopping-energy density (above) and of the chiral 
magnetization (below) for D = 99 (</> model) at the pseudocritical transition point 
J pc (L). 



of the 4 model for {7=1 and D > 0, which emerges from the MC results we shall 
present in this section. We do not know how the two continuous transition lines connect 
to the unique first-order one. There are several possibilities. Under the hypothesis that 
no unique continuous transition occurs, as expected in the IsXY model (see Sec. EJ) , we 
may have: (i) the Ising line becomes of first order at D — < D*, and then the KT 
line meets the first-order line at D* (this possibility is presented in Fig.[HJ); (ii) as before, 
but interchanging Ising and KT lines; (iii) both the Ising and KT lines become of first 
order at D = and D = D sp respectively, with D c ^, D sp < D* , and then they meet at 
D = D*. Our data suggest that the intermediate region between the one characterized 
by two continuous transitions and the one with a unique first-order transition should 
be relatively small, i.e. D ch , D sp , D* are close to 49 and larger than 34. However, 
we cannot exclude that the system undergoes a unique continuous transition at some 
value of D or in some interval of D in this intermediate region. In this case two further 
possibilities may be realized: the Ising and KT line meet and then, a single transition 
line starts, which can be either of first order or of second order (turning eventually to 
first order). 
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Table 6. Results for D = 99 and 49. Here J pc (L) is the pseudocritical J obtained 
by using the equal- weight method, P(-E m in) is the minimum of the normalized energy 
distribution between the two peaks, J m ax(^) is the value of J corresponding to the 
maximum of the specific heat, and C max the maximum of the specific heat. 



L 


J pc {L) 


P(E min ) 


J max (^) 




D = 99 16 


1.5975(3) 


0.120(10) 


1.5956(4) 


0.663(5) 


24 


1.59840(14) 


0.069(23) 


1.59766(13) 


0.600(4) 


32 


1.59942(7) 


0.053(5) 


1.59898(7) 


0.569(2) 


48 


1.59994(4) 


0.032(5) 


1.59973(4) 


0.539(2) 


64 


1.60018(4) 


0.018(3) 


1.60006(4) 


0.522(2) 


D = 49 48 


1.57636(6) 


0.152(13) 


1.57608(6) 


0.283(2) 


64 


1.57661(2) 


0.181(8) 


1.57643(2) 


0.2523(13) 


96 


1.576771(15) 


0.174(14) 


1.576692(15) 


0.2145(13) 


128 


1.576853(11) 


0.151(23) 


1.576801(11) 


0.196(2) 



7.1. First-order transition for large D 

For sufficiently large values of D there is a unique first-order transition from the 
disordered phase to the LT phase. Evidence for a first-order transition is provided 
by the presence of two peaks in the distribution P(E) of the hopping-energy density 
E defined in Eq. (J2DJ). We determine the position of the transition by using the equal- 
weight method for the two phases applied to the energy density. In practice, for each 
lattice size we compute -Emm, the energy corresponding to the minimum between the 
peaks of the histogram. Then, a given configuration is considered as a LT configuration 
if E > Emir,, a HT configuration otherwise. Finally, for each lattice size we determine a 
pseudo-critical value J pc (L) by requiring that for J = J pc (L) the ratio of the probabilities 
of LT and HT configurations is some fixed number w. In our case the value w = 2 seems 
the most appropriate one, since the LT phase is expected to be twofold degenerate. Note 
that J pc (L) converges to the infinite- volume transition value J c for any chosen value w, 
with corrections that are generically of order 1 /V. The position of the transition point 
can also be identified by considering the value J max (£) where the specific heat C takes 
its maximum value. Also J mSLX {L) is expected to converge to J c with corrections of order 
l/V. 

To identify the first-order transition region, we have performed simulations for 
D = 34, 49, and 99. The distributions of the hopping energy and of the chiral 
magnetization for D = 99 shown in Fig. |§] clearly indicate that the transition is of first 
order. They have a deep minimum between the two peaks which decreases exponentially 
as L — > oo, consistently with a first-order transition. For instance, the minimum of the 
energy distribution, P(E min ), is expected to behave as e~ 2aL , where a is the interface 
tension. Fitting the estimates of P{E mm ) reported in Table El to ce~ 2aL , where a is 
the interface tension, we obtain the estimate a = 0.017(2), which is rather small. The 
critical value J c can be estimated by fitting the results for J pc and J max reported in 
Table H to J c + a/V. We obtain J c = 1.60042(5) from J pc and J c = 1.60040(6) from 
<^max- Finally, from the maximum value of C one can also estimate the latent heat 




Figure 11. Hopping-energy distribution for D = 34 and several lattice sizes at the 
corresponding J pc . 



A = E + — E_, where E± are the values of the hopping energy Eh at the peaks of the 
energy distribution. Using A 2 = AC max /V + 0(1/V), we obtain A »s 1.4. 

The analysis of the chiral and spin degrees of freedom does not provide evidence of 
other transitions. On the LT side of the first-order transition, i.e. for J — > J c ~ 1.6004 
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from above, the data for D = 99 reported in Table|H]give rj rs 0.12, which is much below 
the value r] — 1/4 where vortex dissociation destabilizes the LT XY phase. In Table El 
we also report a result for J = 1.6003. This value of J is on the HT side. On the 
other hand, the MC run started from a LT configuration and no drastic change of the 
thermodynamic variables was observed during the run. Thus, the result corresponds to 
an average taken in the metastable phase. Finally, we mention that also the finite-size 
behavior of the derivative of the chiral Binder parameter with respect to J at J pc is 
consistent with the behavior expected at a first-order transition [HI]. Indeed, we find 
dB c (L)/dJ ~ L l l Vcii with l/u cS = 2.0(1), in agreement with 1/ugg = d = 2, which is 
valid in the FSS limit at a first-order transition. 

With decreasing D, the first-order transition becomes weaker and weaker and 
eventually disappears. The results for D = 49 reported in Table El and the distributions 
shown in Fig. El suggest that this occurs for D rs 49. Although two peaks are still 
present, P(E min ) is roughly constant with L and does not show the exponential decrease 
expected at a first-order transition. This is the typical behavior of energy distributions 
close to a tricritical point where the first-order transition line ends; see, e.g., Ref. An 
analogous behavior is observed in the distributions of the chiral magnetization shown in 
Fig- El It is interesting to observe that two peaks in the energy and chiral- magnetization 
distributions are also observed for smaller values of D on lattice sizes of order L ~ 10 2 . 
However, in this case the minimum value of the energy distribution increases. This is 
clearly observed for D = 34, see Fig. ^2 For D = 49 the analysis of the finite-size 
behavior of the derivative of the chiral Binder parameter gives \ jv = 1.84(16). This 
estimate suggests a value smaller than 2, the value expected at a first-order transition, 
and it is consistent with the corresponding value at the tricritical Ising point, which is 
\jv = 9/5 [S3]. We do not find evidence of other transitions for D = 49. 

7.2. Continuous transitions 

7.2.1. Standard finite-size scaling analysis. Let us first discuss the results for D = 1/2. 
In Fig. El we plot R s = £, S /L and R c = £ C /L versus J for L < 500 in the critical 
region. The crossing points correspond to different, though very close values of J: 
Across ~ 1.4688 for R s and J cr0 ss ~ 1.4668 for R c . Apparently, chiral and spin variables 
have different crossing points, as also confirmed by the analysis of the Binder parameters: 
Across ~ 1-4680 for B si p, and J C ross ~ 1.4668 for B c . This difference may be due to scaling 
corrections or to the presence of two different transitions. In order to investigate this 
issue we perform a scaling analysis. We fit the data for 1Z = R s , B^, R c , and B c with 

n(L,J) = f[(J-J c )L 1 /% (42) 

where f(x) is approximated by a polynomial. Fits of the chiral variables R c and B c for 
L > 200 give v = 0.78(2) and v = 0.80(3) with a reasonable, though not good, x 2 ; i- e - 
X 2 /d.o.f. ~ 2. The analysis of the spin variables is much less stable. We obtain v = 0.9- 
1.2 from Bg^ and v ~ 1.5 for R s . These results suggest that spin and chiral variables 
undergo two different, though close transitions. Indeed, if the transition is unique, all 
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Figure 12. ii s = £, S /L (above) and R c = £ C /L (below) versus J for D — 1/2. The 
dotted lines connecting different sets of data are drawn to guide the eye. 



observables have a FSS behavior with the same exponent v. Moreover, the numerical 
results apparently exclude that the chiral transition is in the Ising universality class. 
Indeed, the exponent v appears to be different from the Ising one i/jg = 1 and, at the 
crossing point, R c and B c take the values R c = 1.15(3) and B c = 1.116(5), that differ 
from the Ising values jHHI Rb = 0.9050488292(4), B ls = 1.167923(5). 

Similar results are also obtained for other values of D. For example, the same FSS 
analysis of the data at D = 1/3 up to L = 600 gives J croS s ~ 1-4768 (chiral observables), 
v = 0.79(1) (from R c ), v = 0.80(1) (from B c ), and R c = 1.11(2) and B c = 1.116(5); R s 
and B s ,p have crossing points at J ~ 1.4791 and J m 1.4785, with exponent v = 0.9-1.2. 
For D = 1, R c and B c have a crossing point at J ~ 1.4610, with exponents v = 0.78(2) 
and v = 0.80(2) respectively, and R c = 1.17(3) and B c = 1.132(8); R s and B S( p have 
crossing points at J ~ 1.4632 and J ~ 1.4625, with exponents v = 1.2-1.4. 
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Table 7. Estimates at fixed R c = R\ s of the pseudo-critical coupling Jr c , the chiral 
magnetic susceptibility \c the chiral Binder parameter B c , the derivatives of R c and 
B c with respect to J, the ratio R s and its derivative with respect to J, and the total 
helicity modulus T. They refer to the 4> 4 model with D = 1/2. 



L 


Jr 


AC 


B c 


dR c /dJ 


-dB c /dJ 


R s 


dR s /dJ 


T 


8 


1.5320(9) 


63.76(15) 


1.185(2) 


4 10(71 


0.185(8) 


0.8124(9) 


1.93(3) 


1.470(2) 


12 


1.4976(4) 


107.7(2) 


1.1687(9) 


6.45(6) 


0.340(8) 


0.7798(7) 


2.89(2) 


1.269(2) 


16 


1.4824(3) 


158.8(2) 


1.1659(11) 


9.05(9) 


0.572(10) 


0.7578(7) 


3.94(4) 


1.119(2) 


24 
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32 


1.46766(15) 
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0.7107(7) 


8.28(8) 


0.878(2) 


48 
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13 07(111 
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GO 


1.46499(11) 


1188(3) 


1.1726(8) 


44 8(41 


5.98(10) 


0.6641(7) 


16.87(15) 




G4 
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49.1(5) 


6 87(121 
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696(21 


90 
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1 • 1 KJ 1 JVI J V 7 f 


2323(8) 


1.1811(11) 


76.2(9) 


13.2(3) 


6319(121 

\J *\J\J> 1 ' ■' \ j 


27.5(3) 




120 


1 46523(71 


3790(111 


1 1829(111 


107 4(1 01 

-L KJ I « \ _l ■ V7 y 


21.6(4) 


6080(101 
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1 46528f5l 


4223(91 


1 1855(91 

J . L\J. .7 • 7 \ - 7 / 


117 6(1 01 


24 7(41 


6024(91 
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180 


1.46568(5) 


7622(23) 


1.1910(15) 


185(2) 


45.0(1.0) 


0.5741(12) 


60.2(8) 


0.497(3) 


200 


1.46579(6) 


9132(40) 


1.193(2) 


208(3) 


55.5(1.4) 


0.559(2) 


67.9(1.2) 




256 


1.46601(3) 


14046(37) 


1.1988(13) 


289(3) 


87.3(1.4) 


0.5376(11) 


90.5(9) 


0.421(3) 


300 


1.46624(8) 


18742(148) 


1.202(4) 


364(9) 


119(5) 


0.525(3) 


109(3) 




360 


1.46625(3) 


25534(100) 


1.2002(13) 


438(5) 


159(3) 


0.499(2) 


129.9(1.2) 


0.352(4) 


400 


1.46636(3) 


30907(139) 


1.203(2) 


514(8) 


192(5) 


0.486(2) 


145(2) 




500 


1.46646(3) 


45794(283) 


1.208(3) 


682(13) 


279(9) 


0.462(3) 


183(4) 




512 


1.46651(3) 


48282(338) 


1.201(2) 


696(13) 


283(5) 


0.459(3) 


186(3) 


0.280(6) 


600 


1.46664(6) 


64317(882) 


1.199(3) 


834(30) 


379(18) 


0.439(6) 


209(4) 


0.252(8) 


800 


1.46666(2) 


104255(603) 


1.196(3) 


1095(25) 


564(15) 


0.379(2) 


252(3) 


0.176(5) 


1200 


1.46675(4) 


214664(3746) 


1.186(6) 


1630(91) 


967(57) 


0.302(7) 


317(17) 


0.097(7) 



In conclusion, the analysis presented here apparently favors the presence of two 
transitions. The first one is associated with the chiral degrees of freedom with an 
exponent v ~ 0.8. This is similar to what has been observed in the past in similar 
systems, see Table [D In the next Section, we shall show that the result v ~ 0.8 is a 
crossover effect: The chiral transition belongs to the Ising universality class. 

7.2.2. Finite-size scaling at fixed R c . In order to check whether the appearence of two 
very close transitions is just due to uncontrolled systematic errors, or the value v ~ 0.8 
of the chiral exponent is an effective exponent associated with a preasymptotic regime, 
we performed additional simulations on larger lattices. We follow here the same strategy 
used in Sec. 16. II for the FFXY model, i.e. we perform a FSS analysis at fixed R c = Ri s . 

In Tabled we report the data for D = 1/2 up to L = 1200. The results for R s and 
T at fixed R c , which are also plotted in Fig. 01 and HH1 respectively, strongly favor the 
two-transition picture: R s and T are clearly decreasing and do not show any evidence 

(c) 

of convergence to a nonzero value. Note that the spin correlation length Q at the 
chiral transition is rather large. The results for the largest lattices suggest ^ <: 360, 
explaining the strong crossover effects. The large correlation length ^ also suggests 
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Figure 13. Helicity modulus T at fixed R c for several values of D in the (j) 4 model. For 
comparison, we also show the data for the FFXY model. The dotted lines connecting 
different sets of data are drawn to guide the eye. Note that the data for D = 2 overlap 
with those for the FFXY model. 



that the spin transition should occur at values of J that are very close to J c h- This will 
be confirmed in Sec. 17.2.41 

Similar results are obtained for other values of D, D = 1/20, 1/10, 1/5, 1/3, 1,2,4, 
performing the same FSS analysis for lattices with L < 600. In Fig. El we show 
the helicity modulus T for several values of D. In all cases T appears to decrease, 
consistently with the fact that it should eventually vanish in the large- L limit. 

7.2.3. Ising critical behavior at the chiral transition. We now analyze the critical 
behavior at the chiral transition. In particular, we wish to understand whether the 
results of the FSS analysis of Sec. 17.2.11 indicating a critical behavior different from the 
Ising one, are confirmed by a FSS analysis that includes larger lattices. As in Sec. 16. 1| 
we consider the chiral effective exponents l/v c s obtained from R c and B c . The results 
for D = 1/2 are shown in Fig. EI The curves for l/v e g(L) show first an extended 
region, corresponding approximately to 50 < L < 300, in which w 0.8. Then, the 
FSS curves decrease rapidly towards values which are substantially consistent with the 
Ising exponent v = 1. In Fig. El we show the results for the effective chiral exponent 
r/ cS defined by using Xc ~ L 2 ~ v , in analogy with the definition of u eS . The exponent 
approaches the Ising value rj = 1/4 as L increases. Analogous crossovers are observed 
in other quantities, in the Binder cumulant B c , see Fig. 01 and in the pseudocritical Jr c . 
The Ising regime, in which B c w B\ s + cL _7//4 and Jr c = J c h + aL -11 / 4 , cf. Eqs. ()35j) 
and (|38p. is observed only when L > 600. From linear fits of the data for the largest 
available lattice sizes, we obtain the estimate J c h = 1.4668(1). 

In Figs. El and El we show the chiral Binder cumulant and the effective exponents 
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1/ves derived from R c and B c for other values of D. For < D < 2 the behavior 
is similar to that observed for D = 1/2 and in the FFXY model. In particular, also 
here B c and the exponent v e R are respectively approximately equal to 1.19-1.20 and 
0.8 in a quite large range of values of L. For these values of D we do not observe the 
approach to the Ising value. This is not surprising since, as we shall show in Sec. EJ £s 
increases as D — > and thus our data do not satisfy the condition L ^> ^i c \ For D > 2 
the approach appears qualitatively different and, in particular, no plateau at v w 0.8 
appears. In Table |Hl see also Fig. [THl we give the values of J c h for several values of D. 
In conclusion, a careful FSS analysis on lattice sizes up to L = 1200 shows that the 
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Table 8. Critical values of J at the chiral and KT transitions, as explained in the 



text. 






Jch 




J S p,low from R s 


8 


<P A , D = 


1/3 


1.4765(1) 




1.4792(1) [L = 


600] 


> 0.0018 




1/2 


1.4668(1) 


1.4704(2) 


1.46971(2) [L -- 


= 1024] 


0.0025(2) 


^,D = 


1 


1.4609(1) 




1.4642(1) [L = 


400] 


> 0.0023 


4> A , D = 


2 


1.4696(1) 




1.4736(1) [L = 


400] 


> 0.0027 


4>\D = 


4 


1.4884(1) 




1.4918(1) [L = 


300] 


> 0.0023 


IsXY, C 


= 


1.4684(1) 


1.493(1) 


1.4888(2) [L = 


512] 


0.0167(7) 


FFXY 




2.20632(5) 


2.2415(5) 


2.2348(4) [L = 


1024] 


0.0159(2) 



behavior close to the chiral transition in the 4> 4 model with D = 1/2 is similar to what 
has been observed in the FFXY model. Also in this case the chiral transition should 
belong to the Ising universality class, as expected within the two-transition scenario. 
We expect that, as in the FFXY model, Ising behavior can be clearly observed only 
when L 2> Q ■ Since, d c) Z 360 in this case, much larger lattices would be needed to 
obtain a result as reliable as that of the FFXY model. Similar results are obtained for 
other values of D. In particular, for < D < 2 we observe a rather extended region in 
which u e g w 0.8, as it was already observed in the FFXY model. We will discuss again 
this issue in Sec. El 

7.2.4- The spin transition. In order to determine the critical hopping parameter J sp of 
the spin transition for D = 1/2, we use the same method applied to the FFXY model 
in Sec. 16.21 We assume that the transition belongs to the KT universality class and 
exploit the predicted asymptotic behavior of T and R s = £ s / 1 L at the KT transition, cf. 
Eqs. (jSH| and (J30J). Effective estimates are reported in Table El They allow us to obtain 
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Figure 17. 1/V e ff for several values of D in the </> 4 model, from R c (above) and B c 
(below). 



the final estimate J sp = 1.4704(2). We also mention that the chiral second-moment 
correlation length £c [which can be obtained from the chiral connected correlation 
function (|26|)] is £c = 22(3) at the KT transition. As expected, the Ising and KT 
transitions are very close, indeed 5 = J sp / J c h — 1 = 0.0025(2). 

For other values of D we do not have enough data to perform the previous analysis. 
Nonetheless, we are still able to obtain a lower bound for J sp . We assume only that, in the 
limit L -> oo, R s = 0.7506912... and T = 0.63650817819... for J = J sp , and that R and 
T decrease with L, as suggested by Eqs. (J4T?|) and and confirmed by our data. Then, 
for each L, we determine, using the results for £ S /L and T, the value J for which R s and 
T are equal to the theoretical infinite-volume predictions. These results provide lower 
bounds on the correct critical J sp . If we apply this procedure to the data with D = 1/2 
and L = 1024 we obtain J sPi i ow = 1.46971(2) from £ S /L, and J sp ,iow = 1.46957(2) from 
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Table 9. Critical values of J at the KT transition for the 4> 4 model with D — 1/2. 





L 2 


from R s 


from T 


128 


256 


1.4696(1) 


1.4709(3) 


256 


512 


1.4700(1) 


1.4704(1) 


512 


1024 


1.4705(2) 


1.4704(1) 



T, that are not far from the result obtained above, J sp = 1.4704(2). In TablelHlwe report 
lower bounds for other values of D and estimates of the relative difference 5 = J sp / J c h — 1- 
In all cases S is of the order of 10~ 3 . 

7.3. The phase diagram for small values of D 

Our numerical simulations provide evidence of two transitions for D = 1/3, 1/2, 1,2,4. 
In addition, relatively short MC simulations at D = 1/20, 1/10, 1/5, and 34 appear 
consistent with the two-transition scenario, although the evidence is less robust since 
our simulations used smaller lattice sizes. A natural conjecture is that the two-transition 
scenario extends from D < D*, 34 < D* < 49, down to D — > 0. For D = Hamiltonian 
H,^ is 0(4) symmetric. According to the Mermin- Wagner theorem [HI], it does not 
have any transition at finite temperature, and criticality is observed only for J — ► oo. 
Therefore, the Ising and KT transition lines should meet at D = 0, 1/J = 0, which is a 
multicritical point. 

At a generic multicritical point with two even relevant parameters, such as the 
reduced temperature t = T/T mc — 1 and a generic scaling field g (where T = T mc and 
g = is the position of the multicritical point), the singular part of the free energy is 
given by jOEl EH] 

J- sing ~t 2 - a S(X), X = gt~\ (43) 
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where a and <p are the specific-heat and crossover exponents. This expression is not very 
useful in our case, since T mc = 1/J mc = and £ ~ e cJ . However, it is easy to rewrite 
Eq. (jlBj) in a different form that can be easily generalized. If £(t,g = 0) is the correlation 
length for g = 0, and v is the corresponding critical exponent so that £(t, g — 0) ~ t~ u , 
we can write 



^"sing ~ ^(t, 9 = 0) d A(Y) } Y = g£{t,g = 0)', (44) 

where p = <f>jv. In this expression there is no explicit t dependence and thus we can use 
it for D -> 0. 

In the theory described by 7i^, g is essentially the parameter D. The crossover 
exponent p is related to the RG dimension of the perturbation J2 X ^i x^l x a ^ the 0(4) 
fixed point. In two dimensions, p is equal to the engineering dimension of the operator, 
so that p = 2, with logarithmic corrections. The perturbative analysis is analogous to 
that presented in Ref. [93J, in which the effect of a spin-n perturbation in the 3- vector 
model is discussed. The correct scaling variable is Y = D£ 2 (ln£) -6 , or, using the fact 
that £ ~ J-^V- 7 , 

Y = D(cJ)- 7 exp(2cJ). (45) 

The power appearing in the logarithm, which is universal, has been determined by 
using the perturbative results of Ref. [M] , which gives the anomalous dimension of any 
dimension-zero spin-n perturbation of the N- vector model. In the limit U — > oo, in 
which the <p A model reduces to the 4-vector model with a spin-4 perturbation, we also 
have c = 7r. 

A standard scaling assumption in multicritical theories [HHl ES] is that, along the 
critical lines meeting at the multicritical point, the scaling variables Y and Y are 
asymptotically constant. This implies that the Ising and KT critical lines, corresponding 
to J = Jch(-D) and J = J sp (D), approach the 0(4) multicritical point according to 

-D(cJ c h,sp)" 7 exp(2cJ ch ,sp) = ^ch,sp, (46) 
where Y c ^ and Y sp are two different constants, with Y sp > Y^. It follows 



1 ~ 7 
cJ c h,sp ~ - ln(F chiSp /D) + - In 



~ HY chtSp /D) 



(47) 



and 



i^-l~ h 'W. (48) 
^ch lnY ch /D 

We finally mention that for D < we expect only one KT transition line starting from 
1/J = 0,D = 0, where both fields <f>i tX become critical simultaneously. The same RG 
arguments used above tell us that J sp ,D<o ~ ln(l . 

Unfortunately, we are not able to check numerically these predictions. Indeed, this 
requires to work in a regime in which J c h and J sp are large and this happens only for 
very tiny values of D. But for very small values of D very large lattices are needed in 
order to observe the asymptotic behavior. For D < 1, the values of J c h increase with 
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Figure 19. Sketch of the phase diagram of the Ising-XY model. The continuous, 
dashed, and thick continuous lines represent Ising, KT, and first-order transition lines. 
The distance between the Ising and KT lines is amplified to distinguish them. The 
phase diagram within the two circled regions is unknown. 



decreasing D, see Fig. but D ^ 0.1 seems still far from the asymptotic regime in 
which Eq. (I47J) holds. Note that the presence of the 0(4) point influences the values of 
u c s for the lowest values of D. Indeed, see Fig. El v c r is larger than one for small values 
of L. This is due to the influence of the 0(4) multicritical point, for which 1/u = 0. 

Finally, we would like to point out that an 0(4) zero-temperature multicritical 
behavior is also expected in the antiferromagnetic XXZ model on a triangular lattice. 
Its Hamiltonian is 



where s x is a three-component vector satisfying s x ■ s x = 1. In the easy- axis regime 
A < 1, the ground-state structure is analogous to that of the FFXY model jHZl EE], 
and therefore two close Ising and KT transitions are expected. For A = 1 one recovers 
the isotropic Heisenberg antiferromagnetic model, which is not expected to present any 
finite-temperature transition, but it should become critical only for J — ► oo, see, e.g., 
Refs. |98[ I99j. Due to the effective enlargement of the symmetry from 0(3)® 0(2) to 
0(4) jlOOj . the asymptotic behavior for J —>■ oo is expected to be the same as that of 
the 0(4) vector model. 

8. Phase transitions in the Ising-XY model 

In this section we discuss the critical behavior of the IsXY model (|7|). Some features 
of its phase diagram can be easily determined, see Fig. EH At J = we have an Ising 
transition at C — C\ s = | log(l + y/2) = 0.4406868 . . .. An Ising transition line in the 
J-C plane starts from this point. For C = oo there is a KT transition at jl01[ IHB] 
J pa 1.1199, from which a KT line starts. As we shall see, with decreasing C < C% s the 
KT and Ising transition lines get closer and closer. Another Ising transition occurs at 
J = and C = —C\ s , from which an Ising transition line should start and lie in the 




(49) 
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Table 10. Estimates of several quantities computed at fixed R c = R\ s in the IsXY 
model for C = 0. 



L 


Jr. 


Xc 


B c 


dR c /dJ 


-dB c /dJ 




dR s /dJ 


T 


6 


1.4339(7) 


23.6(2) 


1.1788(6) 


4.74(2) 


1.956(13) 


4.143(4) 


1.973(10) 


0.7886(11) 


8 


1.4415(4) 


38.83(3) 


1.1842(5) 


6.51(2) 


2.822(16) 


5.371(5) 


2.638(10) 


0.7264(9) 


12 


1.4498(3) 


78.02(5) 


1.1904(5) 


10.37(4) 


4.72(3) 


7.680(6) 


3.98(2) 


0.6394(9) 


16 


1.4538(2) 


127.68(12) 


1.1958(7) 


14.73(7) 


6.94(6) 


9.816(14) 


5.42(3) 


0.5780(13) 


24 


1.4589(2) 


257.9(3) 


1.2023(7) 


24.49(12) 


11.93(10) 


13.85(2) 


8.49(5) 


0.4978(16) 


32 


1.46165(7) 


426.7(4) 


1.2043(8) 


35.1(2) 


17.13(14) 


17.58(2) 


11.59(7) 


0.4427(16) 


48 


1.46445(4) 


868.8(1.0) 


1.2061(8) 


58.0(3) 


28.1(2) 


24.25(4) 


17.71(8) 


0.360(2) 


G4 


1.46586(4) 


1441.7(1.6) 


1.2062(7) 


82.1(4) 


39.2(3) 


29.98(5) 


23.38(12) 


0.303(2) 


90 


1.46707(3) 


2627(3) 


1.2029(8) 


121.1(6) 


55.6(5) 


37.55(7) 


31.3(2) 


0.223(2) 


128 


1.46782(3) 


4865(7) 


1.1957(7) 


175.6(1.0) 


76.2(8) 


45.11(10) 


38.9(3) 


0.137(2) 


180 


1.46823(3) 


8855(19) 


1.1843(8) 


240(2) 


97.4(1.2) 


50.29(17) 


41.1(4) 


0.069(2) 


256 


1.46834(3) 


16247(50) 


1.1773(13) 


328(4) 


131(3) 


52.5(3) 


38.7(7) 


0.018(3) 


360 


1.46838(3) 


29466(133) 


1.1724(15) 


463(8) 


177(5) 


52.7(4) 


31.8(9) 


0.005(4) 



region C < — Ci s , bounding an antiferromagnetic phase. We have not investigated the 
behavior of the IsXY model along this antiferromagnetic transition line, although it is 
likely (but we have not numerically checked) that the second-order transition line turns 
into a first-order one as J increases, as it happens for the ferromagnetic case that we 
are going to discuss below. 

The phase diagram of the IsXY model was already studied in Refs. [2E1 EH1 H2j- 
They found an Ising and a KT transition for Ci s > C > 0.1 and a single first-order 
transition for C ^ —4. These works also claimed the existence of unique continuous 
transitions belonging to new universality classes for —4 < C ^ 0.1. Along this line 
critical exponents were found to depend on C. As we shall see, our results do not 
support the existence of a unique continuous transition in which spin and chiral modes 
become critical at the same J. We find two close KT and Ising continuous transitions 
for C ^ —5, and a single first-order transition beyond, in agreement with the argument 
reported in Sec. El which forbids a continuous transition where chiral and spin modes 
become critical. 

We perform the same FSS analysis as described in Sec. 16.11 We first discuss in 
some detail the results for C — 0. In Table ITU1 we report data for C = up to L = 360, 
which turn out to be sufficient to infer the asymptotic critical behavior. As also shown 
in Fig. El the spin correlation length £g at the chiral transition clearly converges to a 
finite value: we find Q = 52.7(4). Correspondingly, the helicity modulus vanishes in 
the large- L limit. Again, the best evidence of an asymptotic Ising critical behavior is 
obtained from the analysis of the chiral Binder cumulant B c . Fig. 0] shows the difference 
AB C = B c — Bjs that clearly converges to the Ising value as L increases. A linear fit to 
a+bL~ 7 ^ of the data for the three largest lattices gives B c = 1.1679(18), to be compared 
with [EEj B lB = 1.167923(5). A linear fit to 6L~ 7/4 of AB C for the three largest lattices 
gives b = 146(6) with x 2 /d.o.f ~ 0.1. The critical coupling Jr c shows the expected rate 
of convergence. A linear fit to a + of the results for the largest lattices gives 
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Figure 20. Effective exponents l/is ft as obtained from the derivatives of R c , B c , and 
R s , for the IsXY model at C = 0. The dotted line corresponds to the Ising value v = 1. 



J c h = 1.4684(1). The effective exponents v e g obtained from the derivative of R c , B c , R s , 
see Fig. EH and the exponent rj ea obtained from Xc behave as observed in the FFXY 
model and in the 4 model in the region < D < 2, see also Sec. El The values of l/^ e ff 
obtained from R c and B c are first larger than the Ising value, then become significantly 
smaller, and eventually converge to 1 from below. 

Estimates of the critical value of J at the KT transition have been obtained following 
the procedure described in Sec. 16 .2\ using the estimate of R s and T with L < 512. We 
obtain J sp = 1.493(1), see also Table El The chiral correlation length ^ at the KT 
transition is £c = 4.3(5). Again, the relative difference 5 = J S p/J c h — 1 of the critical 
couplings is very small, i.e. 5 = 0.0167(7). 

The results for the other values of C are similar. For example, in Figs. |^ E21 and 
ESI we show the helicity modulus T, the chiral Binder parameter B c , and the effective 
exponent l/ves for several values of C. Precise estimates of the spin correlation length 
at the chiral transition are obtained for C = 0.1, 0.2, 0.3: = 26.6(2), 10.8(1), 3.70(2), 
respectively. As expected £s rapidly decreases with increasing C, indicating that the 
difference between the Ising and KT critical temperature becomes larger and larger up 
to C — Cj s . We also performed simulations for negative values of C, C — —0.2, —0.5 
for lattice sizes up to L = 256, C = -1, -2 to L = 128, and C = -3 to L = 192. The 
estimates of R s and T appear to decrease as in the C = case, suggesting that the spin 
modes are not critical at the chiral transition, thus favoring the two-transition picture. 
With decreasing C, ^ apparently increases, but we have not been able to directly 

(c) 

observe that Q is finite on the lattice sizes that we considered. However, our results 
are not consistent with a unique transition. As we shall see in the next section, the 
results for C > — 1 can be explained in terms of an approximately universal crossover 
within the standard two-transition picture. For C = —2 and C = —3 the estimates of 
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Figure 21. Helicity modulus at fixed R c for several values of C in the IsXY model. 
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Figure 22. AB C = B c — B\ s at fixed R c = R\ s for the IsXY at several values of C. 
B\ s = 1.167923(5) is the value of the Binder parameter at the critical point in the Ising 
model ESI. 



RG quantities (see, e.g., Fig. |2Hlib r z^cfr) are apparently stable with L. For instance, the 
analysis of the derivatives of R c , B c , R s with respect to J for C = —2 up to L = 128 
gives \jv « 1.31, 1.40, 1.10. Each estimate appears quite stable with respect to L and 
thus it could be taken as evidence that the asymptotic limit has been reached and that 
the critical chiral behavior is not Ising but belongs to a new universality class. This 
interpretation is however not possible, since, if the transition is unique, the estimate of 
v should be independent of the observable. Analogously, for C = — 3 up to L = 192 we 
obtain \jv w 1.42, 1.57, 1.23 from R c , B c , R s , which are again not consistent with the 
existence of a new universality class. Taking also into account the argument reported 
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in Sec. |2J which forbids continuous transitions where chiral and spin modes become 
critical, we conclude that these results are better explained in terms of crossover effects. 

The FSS behavior of the helicity modulus T at the chiral transition and for lattice 
sizes L <C £g can be related to the fractal dimension of the geometrical Ising clusters. As 
already discussed in Sec. El at the chiral transition spins s effectively live on geometrical 
Ising clusters. If L < &\ XY spins are strongly correlated and thus we can use the 
spin-wave approximation, which leads to the prediction 

T oc L- e (50) 

where e = 2 — d gc > and d gc is the effective fractal dimension of the geometrical clusters 
in this range of lattice sizes. In the case chiral and spin modes are effectively decoupled, 
and therefore the properties of the geometrical clusters are not affected by the spin 
modes, one expects d gc to be equal to the fractal dimension of the geometrical Ising 
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clusters which is predicted to be di gc = 187/96 = 1.9479... [HUES, corresponding to 
e = 5/96 = 0.0521.... For L < the data shown in Fig.ETJshow a power-law behavior 
as in Eq. (JEDJ), but with a rather different exponent. For example, we find e ~ 0.35 for 
0.2 > C > —1. This is due to the interaction between chiral and spin modes, which 

(c) 

changes the effective fractal dimension of the geometrical clusters when L < Q ■ Note 
that the power-law decay of the helicity modulus is another piece of evidence in favor 
of the argument presented in Sec. El which excludes a unique second-order transitions 
for any C, and therefore also for C = —2, —3. 

For larger values of — C, first-order transitions separate the LT from the disordered 
phase. Our simulations indicate the presence of a unique first-order transition for 
C < —5. In particular, we obtain the estimates a ~ 0.012, 0.006 for the interface tension 
at C = —8, —7, respectively. Moreover, for C < —7, the absence of a continuous KT 
transition is rather clear from the analyses of the spin variables. On the LT side of the 
transition, which is at J c = 15.4255(5) for C = —7, we definitely find rj < 1/4: 77 rs 0.035 
at J = 15.426, and rj « 0.043 at J = 15.425 in the metastable LT phase. The presence of 
first-order transitions for C < —5 may also explain the results for the effective exponent 
1/^eff at C — —2, —3. The apparent increase of l/v e s as \C\ increases on small lattices 
can be interpreted as due to the presence of a nearby first-order transition. 

In conclusion, the phase diagram of the IsXY model is characterized by an Ising 
and a KT transition for Ci s > C ^ —5, and by a unique first-order transition for C < C* 
with —5 > C* > — 7. A sketch of the phase diagram is shown in Fig. EIJ Again we do 
not know how the two continuous transition lines connect to the first-order transition 
line. The discussion at the beginning of Sec. [7| applies also in this case, with the only 
difference that a unique continuous transition is excluded by the argument reported in 
Sec. 121 A possible scenario, consistent with our numerical results, is shown in Fig. [T9l 
The crossover curves of the effective exponents are quite similar to those of the </> 4 model, 
at least for C not too negative, as we will further discuss in Sec. H3 

9. Crossover behavior 

In the preceding sections we established that the FFXY model and the 4 and IsXY 
models in a large parameter region undergo two transitions that belong to the Ising and 
KT universality classes respectively. However, the eventual Ising critical regime at the 
chiral transition is reached after a crossover regime which is surprisingly similar in all 
models considered. For example, in most cases we find a rather extended preasymptotic 
region in which the effective correlation-length exponent f e g is approximately equal to 
0.8. Similar values have also been observed in other models related to the FFXY model, 
see Table [D This suggests that this crossover is somewhat universal. 

First, we consider the behavior at the chiral transition and show that any RG 
quantity scales approximately as 

K = ML/l), (51) 
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where / is a model-dependent, but observable-independent length scale. || To 
verify Eq. (j51j) . we consider i? s , determining the factors I so as to obtain 
the best collapse of the data. The results are reported in Fig. |2U The 
data fall on a single curve with remarkable precision. The optimal data 
collapse is obtained for Z/Zffxy ~ 7.0, 4.5, 3.2, 1.8, 1.0, 0.75 respectively for D = 
1/5, 1/3, 1/2, 1, 2, 4, and l/l FFX y = 0.031, 0.089, 0.22, 0.45, 1.11, 2.6, 5.7, 12, 17 for 

|| Note that our results have been obtained at fixed R c and not at the chiral critical point. However, 
since we are dealing with an Ising transition and R c has been fixed to the critical-point Ising value, 
this is irrelevant in the scaling limit. Had we fixed R c to a different value, we would have still observed 
a universal crossover behavior, though the scaling curves would have been quantitatively different. 
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Figure 25. Binder parameters B S( p (top and middle) and B s (bottom) at fixed R c , 
versus L r = L/l, for the FFXY, 4 , and IsXY models. The rescalings I are the same 
as in Fig. El For L r -> oo, B s -> 2 and B s0 -> 3/2. 



C = 0.3, 0.2, 0.1, 0, -0.2, -0.5, -1, -2, -3 in the case of the IsXY model. 

In order to verify the universality of the scaling Ansatz (j5"Tj) we have considered 
the spin Binder parameters B s and B^. In Fig. 122] we plot these quantities in terms 
of L r = L/l, using the rescaling factors determined in the analysis of R s . Again we 
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observe a good universal behavior. % The quality of the data collapse worsens as D or 
— C increases. Indeed, while for R s all data collapse on a single curve, except for a few 
FFXY points that correspond to small values of L (L < 10), the data of B S( j, for D = 4 
(0 4 model) and C = —2, —3 (IsXY model) show some systematic deviations. They 
apparently require a rescaling I that differs by approximately a factor of 2 for D = 4 
and by a factor of 5 for C = —3 from that determined by using R s . Note that these 
deviations are not due to small values of L and/or of Q c , since fi c) « 88 for the 4 
model at D = 4, and > 200 for the IsXY model at C = -2, -3. Most probably, 
the behavior is influenced by the first-order transition line that is present nearby in the 
phase diagram. 

It is easy to realize that, given two different models, I1/I2 corresponds to the ratio 
of the corresponding spin correlation lengths at the chiral transition. Indeed, since 
Rs = fR s (L/l) and £ s — > £s for large L, we must have fn s (x) = a/x for x — > 00, where 
a is model independent, and ^ = l/a. Thus, the scaling we observe implies that R s is 
an approximately universal function of L/^f\ Assuming universality, this allows us to 
estimate the spin correlation length ^ at the chiral transition even in those cases in 

(c) 

which our simulations did not probe the regime L>(] that allows us to obtain a direct 
estimate of ^ ■ For example, using the above-reported estimates of the optimal ratio 
I/Iffxy and = 118(1) for the FFXY model, we infer that « 830, 530, 380, 212, 118 
respectively for D = 1/5,1/3,1/2,1,2 (0 4 model) and £i ~ 131,310,670 respectively 
for C = —0.2, —0.5, —1 (IsXY model). The result for the </> 4 model at D = 2 is consistent 
with the estimate = 116(3) obtained for L = 500. In Fig. 1241 and in those we shall 
present below we have fixed Zffxy = 118 ~ ^i c \ so that we are reporting the data in 
terms of L r = L/^i°\ 

In Fig.|2ni w e do the same analysis for the helicity modulus T. Again, the collapse of 
the data is remarkable. In a quite large region, i.e. for L r < 0.5, T decreases as L~ e with 
e ~ 0.33, with differences of at most 10%. + As already discussed in Sec. El a power- 
law decay is expected in the IsXY model for L ^ and gives the effective fractal 
dimension of the geometrical clusters, d gc = 2 — e. The universal behavior observed in 
Fig. I2IH shows that a power-law decay also characterizes the FFXY and 4 models. It 
is interesting to note that the value e ~ 0.33 is much larger than the value e = 5/96 
corresponding to the fractal dimension d\ gc = 187/96 of Ising geometrical clusters: in 
the IsXY model, this is the dimension of the geometrical clusters only for L ^> ^ s c \ 
This shows, once again, that chiral and spin modes are strongly coupled for L < ^ s c \ A 

% The Binder parameter is expected to converge to 3/2 for L r — > 00 (a simple argument predicts: 
B s <f, = 3/2 + aL r 1 ^ 4 , a > 0). Nonetheless, for L r ps 10, such a quantity is still quite different from its 
asymptotic value. 

+ More precisely, fitting the available data for L > 10 and L r < 0.5 to T = cL~ e , we find e = 0.335(8) 
(x 2 /d.o.f. w 1.1) for the FFXY model, e = 0.332(2) (x 2 /d.o.f. w 0.5) for the <j) 4 model with D = 1/2, 
e = 0.297(4) (x 2 / d -°- f - ~ 0.5) for the </> 4 model with D = 1/5, e = 0.360(5) (x 2 /d.o.f. w 1.1) for the 
IsXY model with C = 0, e = 0.352(3) (x 2 /d.o.f. w 0.7) for the IsXY model with C = -0.5 (d.o.f. is 
the numer of degrees of freedom of the fit). The small differences in the estimates of e can be easily 
accounted for by scaling corrections. 
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Figure 26. Helicity modulus T at fixed R c , versus L r = L/l. The rescalings I are 
the same as in Fig. Note the logarithmic scale on both axes. 



power-law behavior is also observed in R s . As it can be seen from Fig. EU for L r < 0.5 
R s scales as L~ e , although in this case e is much smaller, e ~ 0.1. Note that, if the 
behavior we observe for 0.02 < L r < 0.5 holds up to L r = we would obtain R s , T — > oo 
for L r — > 0. 

In Fig. E?|we plot the effective exponents l/u e s, defined as in Eq. (J37|) . as obtained 
from the derivative of R c , B c , and R s with respect to J, versus the rescaled lattice size 
L r = L/l = L/& ] . We use the same rescaling factors I as above. The agreement is quite 
good, although the quality of the collapse for the chiral exponent z/ e fr is not as impressive 
as what is observed for the spin- related quantities R s , B s , and B s ^. These crossover 
curves show a quite interesting pattern. The effective exponents 1 / u c g derived from R c 
and B c show first a rather flat region in which l/u c s ~ 1.3, corresponding to u e g w 0.8. 
Therefore, a FSS analysis limited to values L < ^ would apparently provide the 
estimate v m 0.8 for the critical exponent associated with the chiral correlation length, 
as found in the FSS analysis of Sec. 17.2.11 Only by using larger lattices, i.e. L 3> ^i c \ can 
one observe the asymptotic Ising critical behavior. Note that, to make things worse, the 
asymptotic approach is nonmonotonic: for L m ^ the exponent u e g starts to increase, 
becomes larger than the Ising value (the curves of l/v e B have a minimum for ~ 2 

corresponding to v c r ~ 1.1 for R c and v e R ~ 1.2 for B c ) and then eventually converges 
to the Ising value. The effective exponents l/u e g associated with the spin variables show 
a similar crossover behavior. The exponent derived from R s must converge to —1 for 
L ^> ^ s c \ since the spin correlation length is finite. However, for L < £s the data of 
1/ueS show a plateau at l/u e s ~ 1.1, see Fig. In the case of B s and B S( p (not shown) 
we observe an analogous plateau, but at l/v e s ~ 1-7. Substantial deviations from the 
above behaviors are observed for D = 4 in the 4 model and for C = —2, —3 in the 
IsXY model, as one can already see from Figs. El and ESI where 1 / u e g is plotted versus 
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Figure 27. Effective exponent l/u e g computed using R c , B c , and R s , versus L r = L/l. 
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L. 

It is impossible to estimate from the data the behavior of the effective exponent u c g 
for L r — > 0, since points are quite scattered. For small values of D in the 4 model l/u e g 
is smaller than one for small lattice sizes, as expected from the influence of the close 
0(4) multicritical point, for which 1/u — 0, see Sec. 17. HI On the other hand, for D > 1 
and C < — 1, the effective exponent approaches the plateau from above, a behavior that 
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Figure 28. Effective exponent -q g at fixed R c , as obtained from \c versus L r = L/l. 
The rescalings I are the same as in Fig. |^J 



may be due to the presence of a first-order transition line for large D and —C. 

In Figs. |2H and EH we plot the effective exponent ?7 e fF obtained from \ c and 
AB C = B c — Bjg. Again, we observe an approximate collapse of the data. A more 
quantitative check can be done by observing that the approximate universality implies 
for L — > oo 

AB C = a(L/£i c) )- 7/4 , (52) 

where a is model independent. If we define b = a(£« ) 7 / 4 , the previous equation gives a 
relation between the ratio of the constants b in different models and the ratio of the spin 
correlation lengths. Using the estimates of b reported in Sees. 16.11 and |H1 for the FFXY 
model and IsXY model with C = 0, we find &ffxy/^i s xy,c=o — 4.0(3). This is in perfect 
agreement with the prediction &ffxy/&i s xy,c=o = (^%fxy/^s%xy,c=o) 7/4 = 4 - 1 ( 1 )- 

The results reported above show that, in the FFXY model and in the <j) A and IsXY 
models for a rather extended region of parameters (respectively for < D < 2 and 
— 1 ^ C ^ 0.2), the finite-size behavior at the chiral transition is somewhat universal. 
A natural conjecture is that this universality is controlled by a multicritical point (or by 
a line of multicritical points of the same type) where chiral and spin degrees of freedom 
become both critical. If this conjecture is correct, one must observe a universal crossover 
behavior also at the spin transition, with the same rescaling factors. 

In Fig.EHwe report R s = C, s /L and T at the spin transitions for the three models for 
which we have a reliable estimate of J sp . The rescaling factors / are those determined at 
the chiral transition. The agreement is again quite good. The existence of scaling at the 
two transitions with the same rescaling factors is another piece of evidence in favor of the 
multicritical origin of the universality we observe. This hypothesis is also supported by 
the fact that the ratio Q / ^ is approximately constant. Indeed, we have ^ = 8.0(5), 
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Figure 29. Binder parameter B c at fixed i? c , versus L r = L/l. The rescalings I are 
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of the Binder parameter at the critical point in the Ising model f£BL 



O c) = 118(1), thus O c) /Cc S) = 15(1), for the FFXY model; ' = 22(3), £f } ~ 380, thus 
£s°V£c^ = 17(3) (assuming an uncertainty of approximately 10% on C,^), for the 4 
model for D = 1/2; £ { c s) = 4.3(5), & ] = 52.7(4), thus & ] /& s) = 12(2), for the IsXY 
model at C — 0. These results are consistent with a constant ratio £s /£c S ^ ~ 15. 

While for L r — ► oo one observes Ising and KT behavior, in the opposite limit 
L r — > 0, the scaling functions converge to the value of the corresponding observable 
at the multicritical point (assuming it exists, of course). The relevant scales are ^ 

(s) 

and £c at the two transitions and thus information on the nature of the multicritical 
point can be obtained by studying the crossover curves in the regimes L <C ffi (chiral 
transition) and L <C Q (spin transition). Unfortunately, in the models in which we 
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transition. On the left we plot the data versus L, on the right versus L r — L/l. The 
rescalings I are the same as in Fig. 1241 



have a precise knowledge of J sp , £c is quite small and thus we are not able to probe the 
multicritical regime at the spin transition. On the other hand, since Q is quite large, 
the behavior at the chiral transition gives us some indications on the multicritical point. 

From the data we can certainly exclude an Ising-XY decoupled multicritical 
point. Indeed, see Fig. |2Hl the helicity modulus is clearly larger than the XY value 
Txy ~ 0.6365. A more likely possibility is that the behavior for L r small is controlled 
by a zero-temperature transition point. Indeed, as we already discussed, our data are 
compatible with R s , T — > oo for L r —> 0. In Fig. EUwe report the results for B s ^ using 
a double logarithmic scale in order to make the region of small L r more visible. Even 
though scaling corrections are present, it is quite plausible that B s< f, — > 1 for L r — > 0, as 
expected in the case of a zero-temperature transition. 



10. Conclusions 



In this paper we investigated the phase diagram of the square-lattice FFXY model and 
of two related models, a lattice discretization of the LGW Hamiltonian for the critical 
modes of the FFXY model, cf. Eq. (£Q), and the IsXY model (JJJ). For this purpose 
we presented a FSS analysis of the numerical results obtained by high-statistics MC 
simulations on square lattices L 2 , for quite large values of L, up to L = O(10 3 ). Our 
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Figure 31. Binder cumulant Bg^ at fixed R c , versus L r = L/l. The rescalings I are 
the same as in Fig. Note the logarithmic scale on both axes. 



main results can be summarized as follows. 

(i) In all models the LT phase is characterized by the breaking of the Z 2 symmetry 
with a chiral magnetization, and by quasi-long-range order associated with the breaking 
of the 0(2) rotational symmetry. We find that the critical behavior of the spin modes 
is controlled by the same line of Gaussian fixed points as in the standard XY model. 
Indeed, we verify that the spin correlation length, the helicity modulus, and the exponent 
r) satisfy the same universal finite-size relations that hold in the LT phase of the XY 
model. 

(ii) We conclusively show that the square-lattice FFXY model undergoes two very 
close transitions, respectively at J sp and J c h, with J sp > J c h- The transition at J sp is 
a KT transition associated with the spin modes, while for J = J c h there is an Ising 
transition where the chiral Z 2 symmetry is restored. For J c h < J < J sp , spins are 
paramagnetic but chiral symmetry is broken. The Ising and KT transitions are very 
close: 5 = (J sp - J ch )/ J ch = 0.0159(2). 

(iii) The phase diagram of the 4> 4 model for U = 1 is shown in Fig. |H1 It presents 
two qualitatively different regions. There is a region in which the behavior is analogous 
to that of the FFXY model, with an Ising and a KT transition at J c h < J sp . These 
transitions are very close: we find 5 = (J sp — J c h)/J c h ~ 0.003 for all values of D in the 
two-transition region (5 = 0.0025(2) for D = 1/2). For D > D* with 34 < D* w 49, 
we find instead a single first-order transition. We do not know how the two continuous 
transition lines connect to the unique first-order one (see the discussion at the beginning 
of Sec.CJ). A possibility is shown in Fig.lHl the Ising line ends at a tricritical Ising point, 
which is the endpoint of the first-order transition line; then the KT line meets the first- 
order one. Since for D = there is no finite-J transition, it is natural to assume that 
the two transition lines get closer and closer as D — > and meet at 1/J = 0,D = 0, 



Multicritical behavior in the fully frustrated XY model and related systems 



50 



which thus acts as an 0(4)-symmetric multicritical point. The critical lines approach 
the multicritical point according to J c h ~ In/} -1 and 5 ~ 1/ InD^ 1 for D — > 0. 

(iv) The phase diagram of the IsXY model (j7J is shown in Fig. EH In Sec. El we 
argue that the IsXY model cannot undergo a continuous transition in which chiral and 
spin modes become both critical. Our MC results support this fact. As in the <j) 4 model, 
there is a rather extented region in which there is an Ising and a KT transition with 
Jch > Jsp- The Ising and KT transitions are close in most of the cases. For example, we 
find 5 = (J sp - J ch )/J ch = 0.0167(7) for C = 0. For C < C*, -7 < C* < -5, we find a 
single first-order transition. Again, we do not know how the two continuous transition 
lines connect to the unique first-order one. A possibility is shown in Fig. EH 

(v) We do not find evidence of continuous transitions — but we cannot completely 
exclude their presence — where chiral and spin modes become both critical, giving rise 
to a new universality class. Thus, our results do not support the field-theoretical results 
of Refs. [701 ESI, which provided some evidence for the existence of a stable fixed point 
in the RG flow of the LGW continuous theory (jJJJ). However, they are not necessarily 
in contradiction, since the lattice models considered here may be outside the attraction 
domain of the stable fixed point found in Refs. [101 ESj- Note also that the existence 
of two continuous transitions at fixed bare parameters (i.e., at fixed D and U, or C) 
means that the critical modes are not appropriately described by the LGW theory, at 
least in the cases considered in this paper. Actually, it suggests that a field theory with 
an Ising and an XY field, for instance the multicritical theory of Refs. [96, 102,, might 
be more appropriate. 

(vi) In most cases the spin correlation length ^ at the chiral transition is quite 
large. For example, ^ ~ 118 in the FFXY model, ^ ~ 53 in the IsXY model for 
C = 0, and ^i c) « 380 in the 4 model for D — 1/2. The presence of this additional 
length scale makes the observation of the expected Ising behavior quite problematic. 

(c) 

We find that Ising behavior is observed only for L > § , which in turn requires very 
large lattices. 

(vii) At the chiral and spin transitions we observe a preasymptotic regime that is 
universal to some extent. Indeed, in the FFXY model, in the <p A model with U — 1 and 
< D < 2, and in the IsXY model with — 1 < C ^ 0.2, we find that RG quantities at 
the transitions scale approximately as 

K(L) n f^L/l), (53) 

where fn(x) is a model-independent function that is of course different at the two 
transitions. The rescalings I are model dependent and set the length scale. Given two 
different models, h/h is equal to the ratios ifi/^fl or £ci/£c2- One might expect 
scaling (j53|l to be valid for L —>■ oo, I — ► oo (or, equivalently, ^i c \^ — ► oo) at fixed 
ratio. However, these conditions are not sufficient to observe this universal behavior. 
Substantial deviations are observed in other parameter regions, e.g., for D = 4 in the </> 4 

(c) 

model and for C = —2, —3 in the IsXY, even though £g is large. This may be due to the 
influence of the first-order transition present for large values of D and —C. Note that, 
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at the chiral transition, the effective chiral exponent v shows an intermediate region in 
which it is approximately equal to 0.8. This explains earlier results obtained by using 

(c) 

smaller lattices. By using only data with L < Q one would obtain v = 0.8 instead 
of the Ising value. Scaling (jSHJ) may be due to a multicritical point where chiral and 
spin modes are both critical. The analysis of the crossover curves for L/l — > allows 
us to exclude a decoupled Ising- XY multicritical transition. A zero-temperature point 
is instead compatible with the numerical data. Note that the universal behavior we 
have observed here concerns the finite-size behavior at criticality, but of course similar 
phenomena are expected in the thermodynamic limit close to the critical points. If 
t c h and t sp are the reduced temperatures and 6 the temperature difference of the two 
transitions, Ising and KT behavior can only be observed if £ c h, t sp <C 5. In the opposite 
limit, 1 ^> t c h,^sp ~ S we expect instead a universal crossover behavior. 

The identification of the multicritical point is a problem that is still open and 
requires further investigation. If the relevant multicritical point is a zero-temperature 
transition, as the MC data suggest, then the 0(4) point discussed in Sec. 17. HI is a 
possible candidate. Otherwise, one should also the consider the finite-temperature 
multicritical point expected in the phase diagram of the generalized IsXY model |2HJ 
with Hamiltonian 

^isXY = - J ^2 ^ A + B(J xO-y)s x ■ s y + Ca x a y ] , (54) 

in frustrated XY models with modulated couplings [TBI E2J D-03 (they are obtained by 
generalizing the FFXY model Q, setting j xy = —r] with rj > on the antiferromagnetic 
links), and in </> theories in which the field-interchange symmetry is broken or in which 
one of the two 0(2) symmetries is broken by adding, for instance, an interaction term 
of the form (<j) 1>x ■ (f> 2 , x ) 2 - 
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Appendix A. Some details on the Monte Carlo algorithms 

Appendix A.l. The FFXY model 

The algorithm we used for the simulation of the FFXY model is based on the results of 
Ref. [HOI- These authors found that the addition of a certain number of overrelaxation 
(OR) sweeps to a local Metropolis sweep makes the algorithm more efficient. They 
also suggested an efficient method to generate the proposed new spin in the Metropolis 



Multicritical behavior in the fully frustrated XY model and related systems 



52 



update. Every sweep over the lattice * one chooses a single unit vector r = (cos a, sin a) 
where a is uniformly distributed in [0, 2tt). Then, for each site x the proposal is 

Sx = 2(f ■ s x )f- s x . (A.T 

This proposal is accepted with the Metropolis probability 



A = Min 



e -{H'-H) 



(A.2) 

The advantage of this proposal compared with more standard ones is that no 
random numbers are needed and also no expensive trigonometric function is evaluated. 
Furthermore, the authors of Ref. [HD] showed that the autocorrelation times in terms 
of sweeps of the Metropolis update with this type of proposal (in the following we will 
refer to it as MR update) are essentially the same as those of a Metropolis update in 
which the new spin is obtained by a random rotation of the old one. At the chiral 
transition the acceptance rate of the MR update is approximately 25%, in agreement 
with Ref. HI]. 

In most of our simulations we used a demon version |lU4j of the Metropolis update. 
We indicate the demon version of any update by a adding a letter D. The demon version 
of the MR update is therefore denoted by DMR. In this framework, one introduces a 
demon field d x 6 [0, oo) with distribution exp(— ^2 x d x ) and replaces the Metropolis 
accept /reject step in the following way. One defines 

d' x = d x -H' + H. (A.3) 

If d' x > 0, the proposal is accepted and s x , d x are replaced by s x , d' x . Otherwise, the 
demon and the spin keep their old value. This way the generation of a random number 
and, more importantly, the evaluation of the exponential is avoided. However, these 
updates do not change the sum 

+ (A.4) 

X 

Therefore, one must also perform heat-bath updates of all demons, by taking d x = 
— In z x , where random number uniformly distributed in (0,1]. It is easy to see 

that by updating d x before the update of s x , one recovers the standard Metropolis 
algorithm. Also note that the acceptance rate of the demon version of the Metropolis 
update is exactly the same as that of the standard one, although the autocorrelation 
times might increase. 

The local overrelaxation (OR) update is given by 

s x = ^ ^ ~ Sx ' (A- 5) 

where S = ^2 y , nnx jxyS y is the sum of s y over the nearest neighbors y of site x. One 
advantage of the OR update is that no random numbers and no evaluations of the 
exponential function are needed. In the case of the standard XY model one finds that 

* We have generated a new r more frequently, namely for each row of the lattice. We have not tested 
whether or not this gives an advantage over choosing a new r just once every sweep. 
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Table Al. Integrated autocorrelation times r c h for the chiral susceptibility of the 
FFXY model close to the chiral transition. r c h is reported in number of measurements. 
# meas. is the total number of measurements and iVi oca i is the number of update cycles 
between two measurements. 



L 


J 


# meas. 


Mocal 


Tch 


256 


2.2053 


22000 


500 


1.76(10) 


384 


2.2062 


56600 


700 


2.67(13) 


512 


2.2063 


39500 


1000 


3.3(3) 


800 


2.2066 


20500 


2000 


3.6(3) 


1000 


2.2063 


35700 


2000 


6.2(5) 



for the optimal mixture of Metropolis and OR updates the dynamical critical exponent 
is reduced to z ~ 1.2 [105J. In contrast, for the FFXY model at the chiral transition 
the authors of Ref. [60] find only a speed-up by a constant factor. 
A complete cycle of the algorithm we used is: 

• Heat-bath of the demons; 

• 5x[(DMR+OR) sweep + iV or OR sweeps]. 

Here (DMR + OR) means that a single-site DMR update is followed by an elementary 
OR update at the same site. At the chiral transition we used N OT = 3, while at the KT 
transition larger values of N OT were used. 

We have implemented a parallel program using the MPI library. To this end, we 
have divided the lattice in one direction into sublattices of size I x L, where I x np = L 
and np is the number of processes that is used. In order to speed up the simulation, we 
also did some kind of parallelization in the second direction, running the OR sweep in 
the second direction in the following way: 

for (xl=0;xl<L/2;xl++) {OR site (xO.xl); OR site (xO,xl+L/2)}, 

where (xO, xl) is a lattice point. This small change allows a better use of the floating- 
point pipelines and gives a speedup of approximately 30%. On a server with 4 1.8GHz 
Opteron CPUs, one complete update cycle takes 0.26s for L = 1000. Roughly speaking, 
one cycle takes 1 s on one Opteron processor with one core. 

In table lAll we report the integrated autocorrelation times of the chiral susceptibility 
for L > 256. The couplings J are close to the chiral critical value J c h = 2.20632(5). 
Fitting all data given in table EU we get N local r ch = 0.022(8) x L L91 ( 6 ) with x 2 /d.o.f.= 
0.99. 

At the KT transition, critical slowing down can be efficiently reduced by using 
overrelaxation updates. In the case of the standard XY model, at the transition one 
should increase N OT as N OT oc L. We found that also in the case of the FFXY model 
at the KT transition this procedure leads to very small autocorrelation times. For 
J = 2.242 [J sp = 2.2415(5)] and L = 1024, by using A" or = 12 we find an integrated 
autocorrelation time r sp pa 5.4(5) in units of cycles (heat bath of the demons followed 
by 5x [(DMR+OR)+12 OR] sweeps). For L pa 1000 the autocorrelation time is smaller 
by a factor of 500 than that at the chiral transition. 
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Table A2. Integrated autocorrelation times r c h for the chiral susceptibility in 
the IsXY model at C = 0, close to the chiral transition. r c h is reported in units 
of measurements. # meas. is the total number of measurements, N\ oca i is the 
number of local update cycles, and A^i ng i e is the number of single-cluster updates 



per measurement. 


L 


J 


# meas. 


Mocal 


single 


Tch 


64 


1.4660 


100000 


500 


20 


0.92(2) 


90 


1.4672 


100000 


700 


20 


1.22(3) 


128 


1.4678 


100000 


1000 


40 


1.62(5) 


180 


1.4681 


100000 


1000 


200 


3.05(12) 


256 


1.4681 


75000 


1000 


100 


6.2(3) 


360 


1.4585 


50000 


2000 


80 


6.1(4) 



Appendix A. 2. The Ising-XY model 

For the IsXY model we used three types of updates. We performed single-cluster jHS] 
and overrelaxation updates of the spins s keeping the chiral field a fixed, and Metropolis 
updates in which a and s are both changed. The proposal for the Metropolis update is 

a' x = — cr x , Sx = (cos a, sin a), (A. 6) 

where a is a random number uniformly distributed in [0, 2tt). Close to the chiral 
transition, the acceptance rate of this update is approximately 11% for C = and 
3% for C = —7. We mostly used the demon version of this update (in the following we 
refer to it as DM update). 

Between two measurements we perform N\ oca x local cycles [a local cycle consists in a 
heat bath of the demons and in 5x(DM sweep + OR sweep)] followed by A^ sing i e single- 
cluster updates. On one 2.2 GHz Opteron CPU, one local cycle takes approximately 
0.20s for L = 512. 

In Table IA2I we report the integrated autocorrelation times of the chiral 
susceptibility for lattices with 64 < L < 360. The couplings J are close to the chiral 
critical value J c h = 1.4684(1). The number of single-cluster updates was chosen ad hoc, 
since the cluster update has little influence on the autocorrelation times. Fitting iV loca ir c h 
for L > 128 (we are neglecting here the role of the single-cluster updates) we obtain 
NiocaiTch = 0.12(4) x L L95 ( 6 ) with x 2 /d-o.f-= 0.53. At the chiral transition, the cluster 
algorithm leaves the dynamic critical exponent substantially unchanged, i.e. z ~ 2. A 
similar observation holds for other values of C, whenever the transition is of second 
order. 

Critical slowing down at the KT transition and in the LT phase is essentially 
eliminated by the single-cluster updates. For instance, for L = 512 and J = 1.492 
we get r sp = 0.94(1) in units of measurements, where a measurement was performed 
after 4 single-cluster updates and 20 local cycles. 
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Appendix A. 3. The 4 model 

Also here, we used a combination of local Metropolis updates and single-cluster updates. 
The single-cluster algorithm was used to update the angle of either 0i or 02. 

We used different versions of the local Metropolis update, depending on the value 
of D. For D < 4 we used a Metropolis update with the following proposal (update M4): 

$H,x = + s ( r ij,x ~ 0.5) , (A.7) 

where (f)j )X = {4>ij, x , 02j,z) and T{^ x are random numbers with uniform distribution in 
[0, 1). The step size s was chosen so that the acceptance rate is roughly 50%. 

This update becomes inefficient for large values of D. Indeed, for D large, 
I 3> | <t>2,x | or vice versa, so that the Metropolis proposal (|A.7|) leads to changes of 
the fields that are of the typical size of the smaller one. Hence, the algorithm becomes 
inefficient. To deal with this problem, we constructed proposals that are specifically 
adapted to this situation: 

(M2) We change only one field at each time. With probability 1/2 we change either 0! 
or 02 using the proposal 

0H )a; = 0iM + s ( r M _o - 5 )> ^'i2,x = 012,3 + s{r 2 , x - 0.5) (A. 8) 

or the analogous one for <p2,x- Here r± tX and r2, x are random numbers with uniform 
distribution in [0,1). The step size s is chosen so that the acceptance rate is 
approximately 25%. This way, we achieve an efficient update of the large component 
of the field. 

(MLE) We exchange the lengths of the two fields. The proposal is given by 

X' l^2,x| ? 7/ I0mI X l A n\ 

I0mI |0mI 

This type of update should speed up the chiral degrees of freedom. The acceptance 
rate is 14% at the chiral transition at D = 49. One could have also exchanged the 
components of the two fields, using 

0M = 0^> 02,x = 0M- ( A -!0) 
However, this proposal has a smaller acceptance rate than (jA.9|) . 

Beside these updates we also used local updates that leave invariant part of the 
Hamiltonian and are close in spirit to overrelaxation updates. A first possibility 
(OR4) consists in performing a field reflection that keeps the 0(4)-invariant part of 
the Hamiltonian constant: 

■>' >lS S-$ x (A.ll) 



S 2 

where S is the sum over the fields of the four neighbors of x. For convenience, we define 
the four-component vector (fi x = (<f>n )X , <f>i2,x, 02M, <t>22,x)- This type of update has no 
free parameters. The remaining part of the Hamiltonian is taken into account in an 
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accept/reject step. We expect the acceptance rate to decrease with increasing D. At 
the chiral transition we find 87%, 70%, 61%, 43%, 40% for D = 1/2,2,4,34,99. Note 
that it seems to have a finite limit as D — > oo. In the simulations we used a demon 
version (DOR4) with a single demon (not one for each site). 

A second possibility is the overrelaxation of one component of the field (OR1). The 
Hamiltonian for a single component <j>ij X of the field, keeping all other fields fixed, has 
the form 

H(<f>ite) = a<j> ijiX + b$ jtX + c<4>- (A. 12) 

The proposal 

-o„., r -<i:!> (A. 13) 

keeps the value of a(f>ij iX + b<f)fj X constant. The quartic part of the Hamitonian is taken 
into account in an accept/reject step. The acceptance rate for this type of move is 67% 
at the chiral transition for D = 49 (56% for D = 0.5). Such a step, in a demon version, 
was performed subsequently for each of the four components of the field at a given site. 
We used a single demon for the whole lattice, refreshing it before the update of each 
site (this is equivalent to using a demon for each lattice site and updating it before each 
sweep). 

Because of the different features of the model for large and small values of D we used 
different local update cycles in the two regimes. For D < 4 we used a local cycle with an 
M4 sweep followed by iV or DOR4 sweeps. In most of the cases we used N OT = 8, but there 
is no sharp dependence of the efficiency on N OT . For large values of D a typical local 
cycle is instead: (M2+DOR4) sweep, 2 DOR4 sweeps, DOR1 sweep, 2 DOR4 sweeps, 
(MLE+DOR4) sweep, and finally 2 DOR4 sweeps. Here (M2+OR4) means a local M2 
update followed by a local DOR4 update on the same site; (MLE+DOR4) has the same 
meaning. Beside the local cycle we also performed single- cluster updates: between two 
measurements we perform Ai ocal local update cycles followed by N siag i e single-cluster 
updates on each field. 

Let us look in more detail at the performance at D = 1/2. At the chiral transition, 
by using the local cycle (M4 + 8 DOR4), we obtain the autocorrelation times reported 
in Table lA"3l They are consistent with a critical exponent z — 2. 

At the KT transition the cluster algorithm allows us to eliminate the slowing down 
of the spin degrees of freedom. For J = 1.47 and L = 2048 we find r sp = 2.2(3) in 
units of measurements. Between two measurements we performed 10 local cycles and 
10 single-cluster updates on each field. 

Appendix B. Some results for the low-temperature phase of 2-d XY models 

In the spin-wave limit the partition function of the 2-d XY model is 

Psw 



z = E/ D ^ exp 



2 

X,fl 
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Table A3. Integrated autocorrelation times r c h for the chiral susceptibility in the 4 
model at D = 1/2, close to the chiral transition. r c h is reported in units of measurents. 
# meas. is the total number of measurements, N\ oca \ is the number of local update 
cycles (one cycle corresponds to M4 + 8 DOR4) and N s i ng \ D the number of single-cluster 
updates for each Geld per measurement. 



L 


J 


# meas. 


Aflocal 


single 


Tch 


360 


1.4663 


100000 


200 


10 


11.0(6) 


512 


1.4665 


43000 


320 


10 


14.5(1.1) 


600 


1.4666 


40000 


200 


10 


40(8) 


800 


1.4665 


127000 


100 


10 


95(11) 


1200 


1.4666 


37600 


200 


10 


>100 



W(n 1 ,n 2 ) J D[0]exp 



SW 



(B.l) 



where n% and n 2 are integer. The shift 2 , Kn fX 8 x l at the boundary takes into account 
winding XY configurations on lattices with periodic boundary conditions. The weights 
are given by 



W(ni, n 2 ) = exp 



-2vr psw [ j^ n i + j^ri 2 



(B.2) 



Note that the contributions from (ni,n 2 ) ^ (0,0) are tiny: for instance, we have 
W(1,0) = 3.487... x 1(T 6 for L x = L 2 at j3 sw = 2/vr. The LT phase of the XY model 
is effectively described by the spin- wave model for fisw > 2 /it. For smaller values of 
Pswi vortices become relevant and disorder the system. Setting s x = cos0 x , the XY 
correlation functions are given by 

(s x s y ) = (B.3) 
J2 ni ,n 2 W{ni,n 2 ) cos[pini(xi - y x ) + p 2 n 2 (x 2 - y 2 )](exrj[i27r((f) x - (f) y )]) G 

where = 2ir/L fl . (.)g denotes the expectation value in a Gaussian system without 
boundary shift. Using the Wick rule 

( {<Px - 4>y) 2n )G = (2n - 1)!! ((4> x - </g 2 >£, (B.4) 
we can express the exponential in terms of the propagator of the Gaussian field: 

(exp^vr^ - (j> y ))) G = exp [-27r 2 ((^ - <P y ) 2 ) G ] . (B.5) 
The exponent rj of the magnetic susceptibility is given by 

The ratio i? = £/L is computed numerically using Eqs. (jB.3|) and (jB.5|) . One can also 
derive the helicity modulus. For a lattice with L\ = L 2 we have 



r = P. 



SW 



^0sw En=-co n " exp(-27T 2 /? w n 2 
J2n=-oo exp(-2ir 2 (3 S wn 2 ) 



(B.7) 
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One may use Eq. (|B.6|) to replace fisw with r\ in Eq. (|B.7|) . and obtain the universal 
relation between the helicity modulus and critical exponent rj in the large- L limit of a 
square L x L lattice. An analogous relation can be obtained for the ratio £/L. 

According to the KT theory, the leading corrections to x> £/L an d T in the LT 
phase are proportional to L~ 2x with x = ir/3 — 2. Using Eq. (jB.6j) . we have 2x = 1/rj — A. 
In the case of \ and £ we also have corrections from the analytic background — they are 
of order L v ~ 2 — which dominate for 2 — rj < 2x. 
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